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1. INTRODUCTION 



The LHC will be a very complex environment with most of the interesting physics signals, and their 
backgrounds, consisting of multi-parton (and lepton/photon) final states. The ATLAS and CMS exper- 
iments will measure these final states with negligible statistical error, even in the early running, and in 
many cases with systematic errors smaller than those achieved by the experiments at the Tevatron (see 
the contribution in these proceedings from G. Dissertori). The luminosity uncertainty and the uncertainty 
in the parton distribution functions (PDFs) can be minimized by the normalization of the physics process 
of interest to certain Standard Model (SM) benchmark processes, such as W, Z, and tt production. Thus, 
it is important to have theoretical predictions at the same or better precision as the experimental measure- 
ments. In many cases, SM backgrounds to non-SM physics can be extrapolated from background-rich to 
signal-rich regions, but a definite determination of the background often requires an accurate knowledge 
of the background cross sections. An accurate knowledge of a cross section requires its calculation to at 
least next-to-leading order (NLO). 

There are many tools for constructing basically any complex final state at the LHC at leading 
order (LO). When interfaced to parton shower Monte Carlo programs, such predictions can provide a 
qualitative prediction of both inclusive and exclusive final states. There are several different interfaces 
between fixed order (both LO and NLO) matrix element and parton shower Monte Carlo programs, with 
a benchmark comparison reported in this workshop. 

A realistic theoretical description of complex final states, though, exists only at NLO 1 , with the 
current limit of such calculations being 2^3 and 2 — ► 4 processes (see below). At LO, calculations 
often have large scale dependence, a sensitivity to kinematic cuts, and a poor modeling of jet structure. 
These deficiencies are most often remedied at NLO. NLO parton level calculations can serve as useful 
benchmarks by themselves, as well as providing an even more complete event description when inter- 
faced with parton shower Monte Carlo programs, or when resummation effects are included. For the 
crucial benchmark processes mentioned above (W, Z and tt production), it is useful to go beyond NLO 
to NNLO. This has been done for W and Z production, including the calculation of differential rapidity 
distributions, and is expected for tt, Z/7+jet and VF+jet production in the near future. Progress towards 
tt is reported in the contributions from M. Czakon, A. Mitov and S. Moch. 

Even at NLO, the calculation of 2 — ► 3 (and 2^4) processes is extremely time- and theorist- 
consuming, so clear priority needs to be established for those processes most needed for the LHC. In the 
2005 Les Houches proceedings, such a realistic NLO wishlist was established (see Table 1). It is grati- 
fying that 3 of the 8 processes (and some which were not listed, for example the one-loop interference 
between gluon fusion and weak boson fusion in Higgs plus dijet production [1,2] ), have been calculated 
in the intervening two years, but daunting to know that 5 remain and a new process has been added. 
As noted in the table, three groups have calculated VFVF+jet since Les Houches 2005 and a detailed 
comparison of the results is presented in these Les Houches proceedings. In addition to the new NLO 
calculation, several processes beyond NLO also have been added to the list. 

The new processes that have been added are: 

• pp — > 6666 

There are several interesting physics signatures involving two b-pairs in the final state, such as 
bbH(—> 66) and hidden valley signatures where Z bosons may decay to multiple 6-quarks. Re- 
lated to this calculation is the production of 4jets, which is less interesting experimentally, but a 
benchmark calculation from a theoretical point of view. 
The calculations beyond NLO added to the 2007 version of the list are: 

• 99 ->■ W*W* 0{a 2 a 3 s ) 

This subprocess is important for understanding the backgrounds for H — ► W^W^*\ 

• NNLO pp -> tt 

'Unless otherwise stated, the terms LO, NLO , NNLO refer to the order in perturbative QCD only. 
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This process is important for the use of tt production at the LHC as a precision benchmark. 

• NNLO to weak vector boson fusion (VBF) and Z/7+jet 

VBF production of a Higgs boson is essential for measuring the coupling of the Higgs to bosons. 
Z/7+jet is an essential experimental process that is used to understand the jet energy scale. It will 
also be useful for PDF determination. 

• In addition, to further reduce the theoretical uncertainty for the benchmark W/Z processes, a 
combined NNLO QCD and NLO electroweak (EW) calculation is needed. The cross sections are 
known separately to NNLO QCD and to NLO EW, but a combined calculation will improve the 
accuracy of the result. 

It is also daunting to realize that all of the three finished calculations from the 2005 list remain 
private code. To be truly useful, such calculations need to be available in programs accessible to exper- 
imenters. Most useful is if the event 4-vectors and event weight outputs can be stored in ROOT n-tuple 
format, so that experimental analysis cuts can be easily applied in a manner similar to what is used for 
the actual data, and so that results do not have to be re-generated if the analysis cuts change. In such a 
format, it is also easy to store not only the nominal event weight, generated with the central PDF of a 
NLO set, but also the weights for the set of error PDFs as well. In such a manner, the PDF uncertainty 
for any event configuration can be easily established, at the expense of a larger n-tuple size. Such a 
modification is being carried out for the MCFM program. 

The calculation of complex multi-parton final states results in the generation of many subtraction 
terms for soft and/or collinear real radiation (e.g. Catani-Seymour dipole or antenna subtraction terms), 
and each of these in turn requires a counter-event to be generated for the Monte Carlo evaluation of the 
matrix element. Thus, for example, in MCFM for W^+2jets (and for Higgs+2jets as well), there are 24 
counter-term events for each real event. The net result is the requirement of a large amount of CPU 
time for computing such cross sections, and the need for many GB of disk space for storing the results 
in ntuples. These requirements will become even more extreme as the complexity of the calculations 
increases. 

Although most of the NLO calculations for multi-particle production so far are private code tai- 
lored to the particular process at hand, there is a clear effort towards more automatisation and making 
results available to the community. Several agreements have been made during the workshop to facilitate 
comparisons and to make at least certain building blocks entering NLO calculations publicly available: 

• Les Houches accord on master integrals: the aim is to have a library of one-loop integrals, finite 
as well as divergent ones, which can be used by anybody using a method which requires scalar 
master integrals. It has been agreed that the format for the labelling of the integrals respectively 
their arguments should follow the LoopTools [13, 14] conventions, as the infrared finite integrals 
are already available in LoopTools. The infrared divergent ones recently have been classified and 
listed in [15] and can be found in analytic form at http : / /qcdloop . f nal . gov. The final 
aim is a webpage containing 

1 . a collection of scalar one-loop integrals in analytic form, 

2. benchmark points and comments which kinematic regions have been tested, 

3. code to calculate the Laurent series of each integral at points specified by the user, 

4. ideally also various codes for the reduction to master integrals. 

This webpage is in Wiki format, such that contributions can be added easily. The location of the 
webpage is 

http://www.ippp.dur.ac.uk/LoopForge/index.php/Main_Page, and input is 
eagerly awaited. 

• If an amplitude is published in an analytic form, numerical values at some benchmark points should 
be given to facilitate cross-checks by other groups. 

All of the 2005 NLO wishlist processes that have been completed to date relied on traditional 
Feynman diagrams for the loop amplitudes. On the other hand, as the complexity of the final-states 
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Process 
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NLO QCD+EW to the VBF channel 
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ZZZ completed by Lazopoulos/ Melnikov/ Petriello [8] 

and WWZ by Hankele/Zeppenfeld [9] 
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relevant for VBF -> H -> VV 
VBF contributions calculated by 
(Bozzi/ ) Jager/ Oleari/Zeppenf eld [10-12] 
various new physics signatures 


NT O ra1 filiations aHHpH to list in 9007 




9. — > 6666 


Higgs and new physics signatures 


Calculations beyond NLO added in 2007 




10. gg 0(a 2 a^) 

11. NNLOpp -> tt 

12. NNLO to VBF and Z/7+jet 


backgrounds to Higgs 

normalization of a benchmark process 

Higgs couplings and SM benchmark 


Calculations including electroweak effects 




13. NNLO QCD+NLO EW for W/Z 


precision calculation of a SM benchmark 



Table 1 : The updated experimenter's wishlist for LHC processes 
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grows further, it may prove necessary to adopt as well new approaches and methods. At the 2007 session 
of Les Houches, several such approaches were under discussion and development, primarily those based 
on the general analytic structure of amplitudes. These methods include recursive techniques at both 
tree and loop level; the use of (generalized) unitarity in four dimensions, and in 4 — 2e dimensions 
(the latter in the context of dimensional regularization); and automated solutions for coefficients of one- 
loop integrals, which is also connected with generalized unitarity. Complex final states possess intricate 
kinematic regions in which either the amplitude itself becomes singular, or a particular representation of 
it becomes numerically unstable. The general identification of such regions, and methods for dealing with 
potential instabilities, are also areas of active interest, which are not unrelated to the use of analyticity to 
construct loop amplitudes. 

Even with the rapid progress we have been seeing in the last few years, there are NLO cross sec- 
tions of interest that will not be completed in a timely manner for the LHC. One question is whether 
we can provide any approximations/estimates of the uncalculated NLO matrix elements based on expe- 
riences with simpler calculations. Table 2 shows the K-factors (NLO/LO) tabulated for some important 
processes at the Tevatron and LHC. Of course, K-factors are a simplified way of presenting the effects 
of NLO corrections (depending on both scale choice and PDF used for example), but the table provides 
some interesting insights. For example, it appears that processes that involve a large color annihilation 
(for example gg — ► Higgs) tend to have large K-factors for scales typically chosen to evaluate the matrix 
elements. The addition of extra legs in the final state tends to result in a smaller K-factor. For example, 
the K-factor for Higgs+2jets is smaller than for Higgs+ljet, which in turn is smaller than that for inclu- 
sive Higgs production. The same is true for the K-factor for VF+2jet being less than that for P^+ljet 
and the K-factor for tt+ljet being less than that for tt. Can we generalize this to estimate that the NLO 
corrections for VF+3jets and tt+2jets will be smaller still? 





Typical scales 


Tevatron i^T-factor 


LHC K-factor 


Process 


Mo 


Ml 


/C(mo) 




/C'(mo) 


K(/io) 






W 


mw 


2mw 


1.33 


1.31 


1.21 


1.15 


1.05 


1.15 


W+ljet 


mw 


Jet 

Pt 


1.42 


1.20 


1.43 


1.21 


1.32 


1.42 


W+2jets 


mw 


jet 

Ft 


1.16 


0.91 


1.29 


0.89 


0.88 


1.10 


WW+jet 


mw 


2mw 


1.19 


1.37 


1.26 


1.33 


1.40 


1.42 


tt 


m t 


2m t 


1.08 


1.31 


1.24 


1.40 


1.59 


1.48 


tt+ljet 


m t 


2m t 


1.13 


1.43 


1.37 


0.97 


1.29 


1.10 


bb 


m b 


2m b 


1.20 


1.21 


2.10 


0.98 


0.84 


2.51 


Higgs 


m H 


jet 

Ft 


2.33 




2.33 


1.72 




2.32 


Higgs via VBF 


m H 


Jet 

Ft 


1.07 


0.97 


1.07 


1.23 


1.34 


1.09 


Higgs+ljet 


m H 


jet 

Ft 


2.02 




2.13 


1.47 




1.90 


Higgs+2jets 


m H 


jet 

P T 








1.15 







Table 2: if-factors for various processes at the Tevatron and the LHC calculated using a selection of input parameters. In all 
cases, the CTEQ6M PDF set is used at NLO. K, uses the CTEQ6L1 set at leading order, whilst K' uses the same set, CTEQ6M, 
as at NLO. For most of the processes listed, jets satisfy the requirements pr > 15 GeV/c and |ry| < 2.5 (5.0) at the Tevatron 
(LHC). For Higgs+l,2jets, a jet cut of 40 GeV/c and \rj\ < 4.5 has been applied. A cut of > 20 GeV/c has been applied 
for the tt+jet process, and a cut of p 3 ^ > 50 GeV/c for WW+jet. In the I4 7 (Higgs)+2jets process the jets are separated by 
AR > 0.52, whilst the VBF calculations are performed for a Higgs boson of mass 120 GeV. In each case the value of the K- 
factor is compared at two often-used scale choices, where the scale indicated is used for both renormalization and factorization 
scales. 
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The dream of experimentalists is for every NLO parton level calculation to come packaged with a 
complete parton shower for the partons produced in the NLO hard scattering process. So far, this exists 
for a few not-too-complicated processes, but it is not so easy to arrange this for each given NLO parton 
level calculation. To make this process easier, it will be useful to have a very systematic shower with a 
simple structure that can be matched to the structure of the NLO calculation. Two programs discussed 
at the workshop, and represented by contributions later in this section, may help. One would naturally 
match to a NLO calculation with antenna subtractions. The other would naturally match to a NLO 
calculation with the widely used Catani-Seymour dipole subtractions. 

For many physics processes, though, we will have to continue to rely upon LO parton shower 
Monte Carlo programs (interfaced with exact LO matrix element calculations). In many instances, a 
large part of the difference between LO and NLO predictions is the use of LO PDFs for the former and 
NLO PDFs for the latter. Nominally, the choice indicated above is correct, but LO PDFs can differ from 
their NLO counterparts by a significant amount due to the influence of DIS data on the global fits. The 
LO PDFs often are changed in such a manner as to lead to significant deviations of LO predictions with 
LO PDFs from NLO predictions with NLO PDFs, in some kinematic regions. One solution that has been 
discussed is the use of NLO PDFs with LO Monte Carlos. This solves the problem mentioned above, 
but can lead to additional problems, for example with predictions for low mass objects at the LHC. The 
solution adopted by several groups, and presented at this workshop, is the development of modified LO 
PDFs, including the best features (for use in LO Monte Carlos) of the LO and NLO PDFs. It will be 
useful/important to tabulate the K-factors using these modified LO PDFs. 

For the maximal exploitation of physics, there are also requirements on the experimental side. We 
suggest that cross sections at the LHC should be quoted at the hadron level, and where possible with 
the estimated parton-to-hadron corrections, so that any theoretical prediction (parton or hadron level) 
can easily be compared after the fact to the archived data [16]. Also, the experimental data needs to be 
quoted only for the range of measurement, rather than extrapolated to the full cross section; for example, 
measurements ofW^eu should be quoted for the range of electron transverse momentum and rapidity 
and of missing transverse energy actually used in the triggering and analysis, rather than performing an 
extrapolation to the full W cross sections. Such recommendations were the exception (CDF M^+jets) 
rather than the rule at the Tevatron and a clear model needs to be set for the LHC. 

The structure of this report is as follows. First a review on expected cross sections and uncertain- 
ties at the LHC from an experimental point of view is given to set the stage. Then various new approaches 
to the calculation of tree-level and one-loop multi-leg amplitudes are presented, followed by a section 
on "improvements on standard techniques", with particular emphasis on the analysis of singularities 
which can create numerical instabilities when integrating multi-particle one-loop amplitudes. Section 
III contains various results, first a tuned comparison of different NLO calculations for pp — > WW+jet, 
then results pointing towards the tt cross section at NNLO, and finally NNLO predictions for hadronic 
event shapes in e + e~ annihilation. The latter is not of direct relevance for the LHC, but is a benchmark 
calculation in what concerns the construction of NNLO Monte Carlo programs in the presence of a com- 
plicated infrared singularity structure. The report is closed by a section on parton showers, addressing 
the matching of parton showers with multi-leg LO matrix elements as well as the matching with partonic 
NLO calculations, which is of primordial interest at present and future TeV colliders. 

2. MEASUREMENTS OF HARD PROCESSES AT THE LHC 2 
2.1 Introduction 

We are approaching the start-up of the world's most powerful particle accelerator ever built. It is expected 
that CERN's Large Hadron Collider (LHC) will start its operation in 2008. Thanks to the unprecedented 
energies and luminosities, it will give particle physicists the possibility to explore the TeV energy range 

Contributed by: G. Dissertori 
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for the first time and hopefully discover new phenomena, which go beyond the so successful Standard 
Model (SM). Among the most prominent new physics scenarios are the appearance of one (or several) 
Higgs bosons, of supersymmetric particles and of signatures for the existence of extra spatial dimensions. 

However, before entering the discovery regime, considerable efforts will be invested in the mea- 
surements of SM processes. We are sure that these have to be seen and thus they can serve as a proof for 
a working detector (a necessary requirement before any claim of discovery is made). Indeed, some of 
the SM processes are also excellent tools to calibrate parts of the detector. However, such measurements 
are also interesting in their own right. We will be able to challenge the SM predictions at unprecedented 
energy and momentum transfer scales, by measuring cross sections and event features for minimum-bias 
events, jet production, W and Z production with their leptonic decays, as well as top quark production. 
This will allow to check the validity of the Monte Carlo generators, both at the highest energy scales 
and at small momentum transfers, such as in models for the omnipresent underlying event. The partem 
distribution functions (pdfs) can be further constrained or measured for the first time in kinematic ranges 
not accessible at HERA. Important tools for pdf studies will be jet+photon production or Drell-Yan pro- 
cesses. Finally, SM processes such as W/Z+jets, multi-jet and top pair production will be important 
backgrounds to a large number of searches for new physics and therefore have to be understood in detail. 

The very early goals to be pursued by the experiments, once the first data are on tape, are three- 
fold : (a) It will be of utmost importance to commission and calibrate the detectors in situ, with physics 
processes as outlined below. The trigger performance has to be understood in as unbiased a manner as 
possible, by analyzing the trigger rates of minimum-bias events, jet events for various thresholds, single 
and di-lepton as well as single and di-photon events, (b) It will be necessary to measure the main SM 
processes and (c) prepare the road for possible discoveries. It is instructive to recall the event statistics 
collected for different types of processes. For an integrated luminosity of 1 fb _1 per experiment, we 
expect about 10 7 W — ► eu events on tape, a factor of ten less Z — > e + e~ and some 10 5 tt — ► n+X events. 
If a trigger bandwidth of about 10% is assumed for QCD jets with transverse momentum pi- > 150 GeV, 
bb — ► fi+X and minimum-bias events, we will write about 10 6 events to tape, for each of these channels. 
Also the existence of supersymmetric particles, for example gluinos with mg « 1 TeV, or a Higgs with 
rriH ~ 130 GeV, would result in sizeable event statistics (10 2 — 10 3 ). This means that the statistical 
uncertainties will be negligible after a few days, for most of the physics cases. The analysis results 
will be dominated by systematic uncertainties, be it the detailed understanding of the detector response, 
theoretical uncertainties or the uncertainty from the luminosity measurements. 

Concerning the experimentally achievable precision, it is worth noting that the numerous quality 
checks during construction and beam tests of series detector modules let us conclude that the detectors as 
built should give a good starting-point performance. Furthermore, cosmic ray muons, beam-gas interac- 
tions and beam halo muons are available as commissioning and calibration tools already before the first 
real proton-proton collisions. Finally, with such first collisions in hand, the trigger and data acquisition 
systems will be timed-in, the data coherence checked, sub-systems synchronized and reconstruction al- 
gorithms debugged and calibrated. The electromagnetic and hadronic calorimeters will be calibrated with 
first physics events. For example, the initial crystal inter-calibration precision of about 4% for the CMS 
ECAL will be improved to about 2% by using the (^-symmetry of the energy deposition in minimum-bias 
and jet events. Later the ultimate precision (« 0.5%) and the absolute calibration will be obtained using 
Z — ► e + e~ decays and the E/p measurements for isolated electrons, such as in W — > eu decays [17]. 
The latter requires a well understood tracking system. The uniformity of the hadronic calorimeters can 
be checked with single pions and jets. In order to obtain the jet energy scale (JES) to a few per-cent 
precision or better, physics processes such as 7 + jet, Z(— ► ££) + jet or W — ► 2 jets in top pair events 
will be analyzed. Finally, the tracker and muon system alignment will be earned out with generic tracks, 
isolated muons or Z — > decays. Regarding all these calibration and alignment efforts, the ultimate 

statistical precision should be achieved very quickly in most cases. Then systematic effects have to be 
faced, which, eg., implies that pushing the tracker R<p alignment from an initial 100 fim to about 10 fim 
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might involve at least one year of data taking. More detailed reviews of the initial detectors and their 
performance can be found in Refs. [18] and [19]. 

The anticipated detector performance leads to the following estimates for the reconstruction pre- 
cision of the most important physics objects : 

• Isolated electrons and photons can be reconstructed with a relative energy resolution characterized 
by a stochastic term (which is proportional to 1/y/E) of a few per-cent and an aimed-for 0.5% 
constant term. Typically isolation requirements are defined by putting a cone around the elec- 
tron/photon and counting the additional electromagnetic and hadronic energy and/or track trans- 
verse momentum within this cone. The optimal cone size in rj — 4> space 3 depends on the particular 
analysis and event topology. For typical acceptance cuts, such as a transverse momentum above 
10-20 GeV and \rj\ < 2.5, electrons and photons can be expected to be reconstructed with ex- 
cellent angular resolution, high efficiency (> 90%) and small backgrounds. Again, the precise 
values depend very much on the final state topology and the corresponding tightness of the selec- 
tion cuts. Most importantly, the systematic uncertainty on the reconstruction efficiency should be 
controllable at the 1-2% level, using in-situ measurements such Z — ► e + e" decays, with one of 
the electrons serving as tag lepton and the other one as probe object for which the efficiency is 
determined. 

• Isolated muons, with similar acceptance cuts as mentioned above for electrons, should be recon- 
structed with a relative transverse momentum resolution of 1 - 5% and excellent angular resolution 
up to several hundreds of GeV. Again, a systematic uncertainty on the reconstruction efficiency of 
1-2% appears to be achievable. 

• Hadronic jets will be reconstructed up to pseudo-rapidities of 4.5 - 5, with good angular resolu- 
tion. The energy resolution depends rather strongly on the specific calorimeter performance. For 
example, in the case of ATLAS (CMS) a stochastic term of the order of 50 - 60% (100 - 150%) 
is to be expected when energy deposits in projective calorimeter towers are used for the jet clus- 
tering procedure. Important improvements on the CMS jet energy resolution are expected from 
new approaches such as particle flow algorithms. Well above the trigger thresholds jets will be re- 
constructed with very high efficiency; the challenge is the understanding of the efficiency turn-on 
curves. In contrast to leptons, for jets the experimental systematic uncertainties are much more 
sizeable and difficult to control. A more detailed discussion will follow below. 

A further important question is the lowest px threshold above which jets can be reconstructed 
reliably. Contrary to the naive expectation that only high-px objects (around 100 GeV and higher) 
are relevant, it turns out that many physics channels require jets to be reconstructed with rather 
low transverse momentum of ~ 20 — 30 GeV. One reason for this is the importance of jet veto 
requirements in searches for new physics, such as in the H — ► WW* — > 21 2u channel, where 
a jet veto is necessary to reduce the top background. The experimental difficulties related to the 
understanding of the low -pr jet response 4 , the thresholds due to noise suppression, the impact of 
the underlying event and additional pile-up events and ultimately the knowledge of the JES lead to 
the conclusion that it will be extremely challenging, if not impossible, to reliably reconstruct jets 
below a pr of 30 GeV. In addition, also the theoretical predictions are challenged by very low-px 
effects, as for example induced by jet veto requirements. Here fixed-order calculations may have 
to be supplemented by resummations of large logarithms. 

• Finally, the missing transverse energy will be a very important "indirect" observable, which is 
constructed from measurements of other quantities, such as all calorimeter energy deposits. Many 
searches for new physics, such as Supersymmetry, rely very much on this observable. However, 
it turns out that it is also an extremely difficult quantity to measure, since it is sensitive to almost 
every detail of the detector performance. Here it is even more difficult to give estimates of the 

3 Here 77 denotes the pseudo-rapidity and <f> the azimuthal angle around the beam pipe. 
4 The jet response is defined as the ratio of the reconstructed and the "true" jet momentum. 
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expected systematic uncertainties. Also, the reconstruction performance depends very much on 
the details of the particular final state, such as the number of jets and/or leptons in the event, the 
existence of "true" missing energy, e.g. from neutrinos, the amount of pile-up events and in general 
the overall transverse energy deposited in the detector. The very first data will be of paramount 
importance for a timely understanding of this quantity. 

More detailed discussions of the expected detector and reconstruction performance can be found in recent 
reviews ( [18], [19]), for ATLAS in Ref. [20] and for CMS in its Physics Technical Design Reports 
(PTDR), Vol. 1 [17] and Vol. 2 [21]. 

In the following I will concentrate on the early physics reach of the LHC experiments, i.e. on 
measurements to be performed on the first few hundred pb" 1 up to 1 fb _1 of integrated luminosity. Many 
reviews exist on this topic, such as Refs. [19,22-24] to mention only a few. Most of the results presented 
here are taken from the CMS PTDR Vol. 2 [21], because it represents the most recent comprehensive 
overview compiled by one of the LHC experiments. 

2.2 Jet production 

Because of its extremely large cross section, the inclusive dijet production (pp — ► 2 jets + anything) 
completely dominates over all other expected LHC processes with large momentum transfer. At low- 
est order in perturbative Quantum Chromodynamics (QCD), it is described as a 2 — ► 2 scattering of 
partons (quarks and gluons), with only partons in the initial, intermediate and final state. Depending 
on the exchanged transverse momentum (or generally the energy scale of the scattering process), the 
final state will consist of more or less energetic "jets" which arise from the fragmentation of the outgo- 
ing partons. Indeed, soft scattering processes, which give the largest contribution to the total inelastic 
proton-proton cross section, are most likely, leading to final states with hundreds of soft (i.e. below a few 
GeV) charged and neutral hadrons, uniformly distributed over most of the experimental acceptance in 
pseudo-rapidity. Since these are the most likely processes to occur, they are triggered on with the least 
stringent requirements and thus called "minimum-bias" events. For the same reason they also represent 
the typical pile-up events which can occur simultaneously with other triggered proton-proton collisions. 
Therefore very early measurements of the production rates 5 and the charged particle distributions will be 
extremely important, in particular for the tuning of the widely used Monte Carlo generators. Here I will 
not discuss further this class of measurements, but rather concentrate on the parton scattering at large 
transverse momentum. Examples of envisaged studies of minimum-bias events can be found in [21]. 

For outgoing partons with transverse momentum well above the QCD fragmentation scale (A ~ 
1 GeV) the picture of jet production arises, namely well collimated bundles of particles, leading to 
isolated clusters of deposited energy in the calorimeters. Several algorithms exist for the clustering 
of the final state objects (simulated particles, calorimeter towers, charged tracks) into jets with a well 
defined four-momentum, which in the optimal case closely matches the four-momentum of the original 
scattered parton. Examples of commonly used prescriptions are the Iterative Cone, Midpoint Cone, 
SISCone and kx algorithms. In particular, the latter two algorithms recently receive a lot of theoretical 
and experimental attention, mainly because of their property of being infrared and collinear safe to all 
orders of perturbation theory. A detailed discussion of those jet algorithms is given elsewhere in these 
proceedings, as well as in [25-27] and references therein. 

For the measurement of the inclusive jet cross section we simply count the number of jets inside 
a fixed pseudo-rapidity region as a function of jet pr- For a second typical measurement, the dijet cross 
section, events are selected in which the two highest px jets, the leading jets, are both inside a specified 
pseudo-rapidity region and counted as a function of the dijet (invariant) mass. Both cases are inclusive 
processes dominated by the 2^2 QCD scattering of partons. The distinction between inclusive jets and 

5 Currently the extrapolations from the TEVATRON up to the LHC energies suffer from large uncertainties. For example, 
various Monte Carlo generators predict charged track multiplicities which differ by more than 30%. 
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dijets is only in a different way of measuring the same process. For a common choice of the i] region, 
events selected by the dijet analysis are a subset of the events selected by the inclusive jet analysis, but 
the number of events in the two analyses coming from QCD is expected to be close at high px- The 
steeply falling cross sections are shown in Fig. 1. For the inclusive jet case, the spectrum roughly follows 
a power law, however, with increasing power for increasing px, ie., the power increases from about 6 at 
Pt = 150 GeV to about 13 at px = 3 TeV and keeps on increasing with jet px- 
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Fig. 1: Inclusive jet (left) and dijet (right) cross section measurements as foreseen by CMS [21]. The central cross section 
values are taken from a leading-order calculation in dependence of the transverse momenta of the hard interaction. The insert 
on the right plot indicates various trigger paths. 



It can be seen that even for very small integrated luminosities the statistical uncertainties will be 
negligible, up to very high jet momenta. Thus the TEVATRON reach in terms of highest momenta and 
therefore sensitivity to new physics, such as contact interactions or heavy resonances, will be quickly 
surpassed. For 1 fb^ 1 , the inclusive cross section for central jet production (ie. jet pseudo-rapidities 
below ~ 1) will be known statistically to better than 1% up to a px of 1 TeV, and the statistical errors on 
the dijet cross section will be below 5% up to dijet masses of 3 TeV. 

The real challenge for these measurements will be the determination and control of the jet energy 
scale. As mentioned above, the cross sections are steeply falling as a function of jet px- Therefore any 
relative uncertainty on the jet px will translate into a n-times larger relative uncertainty on the cross 
section, where n indicates the power of the spectrum in a specified px region, ie. da/dpx oc p T n . For 
example, a 5% uncertainty on the energy scale for jets around 100-200 GeV of transverse momentum 
induces a 30% uncertainty on the inclusive jet cross section. This is also shown in Fig. 2 (left), here 
for the case of a 3% JES uncertainty. As a comparison, in Fig. 2 (right) we see the expected theoretical 
uncertainties on the inclusive jet cross section from the propagation of pdf uncertainties. These are below 
the 10% level up to a jet px of 1 TeV, thus much smaller than the experimental systematics from the JES. 
Therefore it is obvious that a measurement of the inclusive jet cross section will not allow to constrain the 
pdfs, unless the JES is known to 2% or better. This is definitely beyond reach for the early phase of the 
LHC, and might remain a huge challenge even later. Furthermore, because of these large experimental 
uncertainties, it might turn out that the currently known next-to-leading order (NLO) perturbative QCD 
calculation of the hard scattering process is precise enough for a comparison to data. However, with 
better experimental control at a later stage and/or other definitions of observables (see below) the need 
for going to next-to-next-to-leading order (NNLO) might arise. 

Obviously, the knowledge of the JES also has a strong impact on the achievable precision of the 
dijet cross section measurement, as shown in Fig. 3 (left). However, the problem can be avoided by 
performing relative instead of absolute cross section measurements. A well suited observable is the dijet 
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Fig. 2: Left: Relative systematic uncertainties of the inclusive jet cross sections for the kr algorithm versus jet pr due to 
a change in the JES of ±3% for three bins in rapidity, y. The error bars indicate the statistical uncertainty. Right: relative 
uncertainties propagated from the error sets of the CTEQ6M [28] pdfs, for the same regions in rapidity. Plots taken from [21]. 



ratio N(\rj\ < |??in|)/-^V(|r/i n | < \t]\ < |??out|), ie., the ratio of the number of dijet events within an 
inner region \r]\ < \r]- m \ to the number of dijet events within an outer region \r/- m \ < \rj\ < \r] out \. Both 
leading jets of the dijet event must satisfy the \?]\ cuts. In Ref. [21] the values chosen were r] m = 0.5 
and ?] out = 1, whereas in a recent update [29] of the CMS studies on inclusive and dijet production 
they have been increased to 0.7 and 1.3, respectively. The dijet ratio has two interesting features. First, 
it is very sensitive to new physics, such as contact interactions or the production of a heavy resonance, 
because those lead to jets at more central rapidities than in genuine QCD dijet events. Second, in the 
ratio we can expect many systematic uncertainties to cancel. For example, the luminosity uncertainty 
completely disappears in the ratio. More importantly, also the JES uncertainty is strongly reduced, since 
the dijet ratio is sensitive only to the relative knowledge of the scale as a function of rapidity, but not to 
the absolute scale any more. This is well illustrated in Fig. 3 (right), where the JES uncertainty is shown 
to be reduced to about 3%. In this figure also the sensitivity to new contact interactions at various scales 
is indicated. Hence we have a nice example of a ratio measurement where systematic uncertainties are 
reduced. Having an observable in hand with experimental systematic uncertainties at the level of 5% or 
less, it might become relevant to obtain a NNLO prediction for jet production. 

As we have seen above, the JES is the dominant source of uncertainty in jet cross section measure- 
ments. Obviously, it is also important for many other analyses and searches which involve jet final states 
and possibly invariant mass reconstructions with jets. Therefore major efforts are devoted by the ex- 
perimental collaborations to prepare the tools for obtaining JES corrections, both from the Monte Carlo 
simulations and, more importantly, from the data themselves. Currently approaches are followed which 
are inspired by the TEVATRON experience [30, 31]. The correction procedure is split into several steps, 
such as offset corrections (noise, thresholds, pile-up), relative corrections as a function of i], absolute 
corrections within a restricted ^-region, corrections to the parton level, flavour-specific corrections etc. 
At the LHC startup we will have to rely on Monte Carlo corrections only, but with the first data coming 
in it will be possible to switch to data-driven corrections. At a later stage, after a lot of effort will have 
gone into the careful tuning of the Monte Carlo simulations, it might be feasible to use Monte Carlo 
corrections again. A rough estimate for the early JES uncertainty evolution in CMS is 10% at start-up, 
7% after 100 pb _1 and 5% after 1 fb _1 [32]. Certainly it will be difficult and require time to obtain a 
detailed understanding of the non-Gaussian tails in the jet energy resolution. 
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Fig. 3: Left : Systematic uncertainty on the dijet cross section due to JES (solid curve), parton distributions (dashed curve) and 
calorimeter energy and r) resolution (dotted curve), compared to the statistical uncertainties for 10 fb -1 (errorbars). Right : 
Systematic bounds on the dijet ratio from uncertainties in the relative JES (dashed curve), parton distributions (dotted curve) 
and calorimeter energy and r\ resolution (dot dashed curve), compared to the expectations of QCD and three contact interaction 
scales (solid line and curves). Plots taken from [21]. 



Concerning data-driven JES corrections, one of the best channels is 7+jet production. At leading 
order, the photon and the jet are produced back-to-back, thus the precisely measured photon energy can 
be used to balance the jet energy. Real life is more difficult, mainly because of additional QCD radiation 
and the large background from jets faking a photon. These can be suppressed very strongly with tight 
selection and isolation cuts (eg., no additional third jet with a transverse energy beyond a certain threshold 
and tight requirements on additional charged and neutral energy in a cone around the photon). The need 
to understand well the photon-faking jet background and the photon fragmentation is avoided by using 
the channel Z(— > ^)+jet, with electrons or muons, however, at the price of a lower cross section. 

Besides being a tool for obtaining JES corrections, both 7+jet and Z + jet processes will also be 
important handles for constraining the gluon pdf. It appears feasible to probe the gluon pdf at Bjorken-x 
values between about 0.0005 and 0.2 with a few per-cent statistical errors after only 1 fb _1 of integrated 
luminosity [33]. The x value is well determined using the lepton or photon kinematics only, thus it does 
not suffer from the less precise measurement of the jet momentum. Of course, in order to consistently 
constrain NNLO pdf sets (which should become more and more relevant with time), a NNLO calculation 
of the hard scattering part of the process is needed. Whereas this appears beyond reach for the 7+jet 
case, the Z+jet process might be tractable within the not-too-far future. As discussed below, Z+jet (as 
well as W+jet) production is a very important background to many searches, therefore having a NNLO 
prediction should be very valuable, also as a benchmark for Monte Carlo generators which combine 
leading order (LO) and/or NLO matrix elements with parton shower models. 

2.3 Vector boson production 

The production of vector bosons (W and Z), triggered on with their subsequent leptonic decays, will be 
among the most important and most precise tests of the SM at the LHC. The leptonic channels, mainly 
electrons and muons, can be reconstructed very cleanly, at high statistics, with excellent resolution and 
efficiency and very small backgrounds. At the same time, the theoretical predictions are known to high 
accuracy, as discussed in more detail below. This precision will be useful for constraining pdfs, by 
measuring the rapidity dependence of the Z production cross section, in particular when going to large 
rapidities and thus probing low x values. As proposed in [34], this process will serve as a standard candle 
for determining to high precision (at the few per-cent level) the proton-proton luminosity or alternatively 
the parton-parton luminosity. Finally, it will be attempted to improve on the current precision of the W 
mass. Besides that, W and Z production will be an important experimenter's tool. As mentioned already 
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earlier, Z and W decays to leptons will be used to understand and calibrate various sub-detectors, measure 
the lepton reconstruction efficiencies and control even the missing transverse energy measurement. 

Below I will first discuss the inclusive case, concentrating on resonant production. Then I will 
highlight some issues for the W and Z production in association with jets. Although being highly in- 
teresting processes, di-boson production will not be discussed here, since for integrated luminosities up 
to 1 fb" 1 the statistical precision will be the limiting factor for these measurements and only allow first 
proofs of existence and rough validations of the model expectations. 

2.3.1 Inclusive W and Z production 

Inclusive W and Z production currently is and probably will remain the theoretically best known process 
at the LHC. Predictions are available at NNLO in perturbative QCD, fully differential in the vector 
boson and even the lepton momenta [35]. Figure 4 (left) shows the Z rapidity distribution at various 
orders in perturbation theory. We see that the shape stabilizes when going to higher orders and that the 
NNLO prediction nicely falls within the uncertainty band of the NLO expansion, giving confidence in 
the good convergence of the perturbation series. More importantly, the renormalization scale uncertainty 
is strongly reduced at NNLO, to a level of about 1% for Z rapidities below 3. A renormalization scale 
uncertainty even below 1% can be obtained for ratio observables such as a(W + )/<7(W~ ) and <r(W)/cr(Z), 
possibly as a function of rapidity. Again, ratio measurements are interesting also from the experimental 
point of few, since many systematic uncertainties cancel completely or to a large extend. The prospect 
of a precise measurement and knowing the hard scattering part of the process so well means that we 
have a tool for precisely constraining pdfs (or couplings and masses, in a more general sense). Indeed, 
when taking the full theoretical prediction for the W and Z production cross section, ie., the convolution 
of pdfs and hard scattering part, its uncertainty is dominated by the limited knowledge of the pdfs, 
currently estimated to be around 5-7% [36, 37]. This will then also limit the proton-proton luminosity 
to a precision of this size, unless the pdfs are further constrained, mainly by the rapidity dependence of 
the cross section, as for example shown in Ref. [36]. It is worth noting that at this level of precision also 
electro-weak corrections have to be considered [38-40]. 

An important point to make in this context is the importance of having differential cross section 
predictions. If we take resonant W and Z production at central vector boson rapidity, we probe x values 
of around 0.006, a region rather well constrained by the current pdf fits. However, for larger rapidities 
we probe more and more the small x region, which is less well known, eg., at leading order and for a 
Z rapidity of 3 we need (anti-)quark pdfs at x = 0.12 and x = 0.0003. Experimentally, because of 
the detector acceptance, we can only access a limited sub-region of the full phase space. This means 
that when measuring a total cross section, we have to extrapolate the measurement to the full accep- 
tance (eg., full rapidity), which introduces a model dependence, especially on the poorly known low-x 
region. On the other hand, having differential predictions, we can compute exactly the same quantity as 
we measure, thus eliminating any extrapolation uncertainty. Similarly, for constraining NLO (NNLO) 
pdfs, exactly the same acceptance cuts (on the leptons) as in the data can now be applied on the avail- 
able NLO (NNLO) predictions. Of course, with more and more differential higher-order predictions 
becoming available, this kind of argument applies to any cross section measurement (and/or deduced 
determination of physics quantities such as couplings, masses, pdfs), namely that we should compare 
the measurements and predictions for the experimentally accessible acceptance and avoid un-necessary 
extrapolations, which will not teach us anything new and only introduce additional uncertainties. 

As mentioned above, the experimental reconstruction of W and Z production is rather straight 
forward. Leptons are required to have a minimum pt of about 20 GeV, within a pseudo-rapidity of 2.5 
(cf. Fig. 4, right). In the Z case the mass peak allows for further event selections and background estima- 
tions. However, the neutrino in the W decay leads to missing energy, which obviously is reconstructed 
less precisely. Instead of an invariant mass peak only the transverse W mass can be reconstructed, with 
larger backgrounds than for the Z. Here it is interesting to mention that a jet veto can help to control bet- 



14 



> 

O 



ft 



a 

13 



b 

CM 




24000 
22000 
20000 

03 

18000 
"£16000 
■Z 14000 
° 12000 
-Q 10000 
§ 8000 
z 6000 
4000 
2000 
0, 




generated Z 

generated Z with (nfl^^ . 
generated Z (bb) 
reconstructed Z (bb) 
reconstructed Z (be) 
reconstructed Z (ee) 



o 1 



3 4 5 



8 9 



Fig. 4: Left : QCD predictions at various orders of perturbation theory for the Z rapidity distribution at the LHC. The shaded 
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ter the QCD backgrounds and to improve the resolution of the missing transverse energy reconstruction. 
However, a jet veto introduces sensitivity to low-p^ QCD radiation, thus comparing the measurement to 
a calculation for the same acceptance cuts will only be meaningful if soft-p^ resummation effects are 
taken into account in the predictions. Fortunately, with the Z+jet process we have an experimental han- 
dle to study these issues rather precisely (see also below), since the radiation pattern in W+jet and Z+jet 
events is very similar. In Ref. [21] it has been shown that reconstruction efficiencies and ultimately cross 
section measurements with systematic uncertainties around 2% (or better) should be possible, excluding 
the luminosity uncertainty. 

2.3.2 W/Z+jets production 

Vector bosons produced in association with jets lead to final states with high-p^ leptons, jets and possibly 
missing transverse energy. Such a topology is also expected for many searches, in particular for squark 
and gluino production and subsequent cascade decays. Obviously it will be important to understand 
these SM processes as quickly as possible and validate the available Monte Carlo generators, which typ- 
ically combine LO matrix elements with parton showers. A standard observable will be the W/Z cross 
section as a function of the associated leading jet transverse momentum or the number of additional jets. 
Obviously, such measurements will suffer from the same JES uncertainties as the QCD measurements 
discussed above, and thus constitute only limited calibration tools during the early data taking. The prob- 
lem can be reduced by defining clever ratios of cross sections, involving different vector bosons and/or 
number of additional jets, or by normalizing the predictions to the data in limited regions of the phase 
space (eg. for small jet multiplicity and extrapolating to larger multiplicities). A completely different 
approach is to take a more inclusive look at this process, in the sense that the Z transverse momen- 
tum is measured from the lepton kinematics, which is possible at high statistical and, more importantly, 
high experimental accuracy (cf. Fig. 5). This distribution can be understood as the convolution of the 
Z+0/1/2/. . .jets distributions, therefore any model intended to describe Z+jets production has necessarily 
to reproduce the Z px distribution over its full range. As mentioned above, in this context it would be 
highly desirable to have a NNLO prediction, possibly matched with a resummation calculation, for a 
comparison to the precise data and as benchmark for other approximations, implemented in Monte Carlo 
simulations. 



15 




Fig. 5: Expected transverse momentum distribution for Z production at the LHC, for a pr range up to 20 GeV (left) and 100 
GeV (right). The pr of the Z is reconstructed from the lepton kinematics. The fluctuations in the spectra indicate the statistical 
precision achievable with 1 fb _1 of integrated luminosity. Plots taken from [36]. 



2.4 Top pair production 

The top quark is produced very abundantly at the LHC. With 1 fb _1 of integrated luminosity, we should 
already have a couple of thousand clean signal events on tape in the di-lepton channel, and a factor of 
10 more in the single lepton channel (lepton+jets channel) [21]. The physics case for the study of top 
production is very rich and can not be discussed in detail here. For example, a recent review can be found 
in Ref. [42]. Combining many different channels, a top mass measurement with a precision of 1 GeV 
might be achieved, which together with a precise W mass measurement constitutes an important indirect 
constraint of SM predictions and its extensions. The production cross section (for single and top-pair 
production) will be an important measurement, again for testing the SM predictions and because top 
production is a copious background to a large number of new physics searches. In the single muon+X 
channel, the top-pair production cross section will soon (ie. with about 1 fb _1 ) be measured with a 
statistical precision of 1%. The total uncertainty of 10-15% (excluding the luminosity uncertainty) will 
be dominated by systematics, most notably due to the knowledge of the b-tagging efficiency. At the 
moment it seems difficult to reduce this uncertainty to below 10% [21], even for much larger integrated 
luminosities. Therefore this should be seen as a benchmark value to be challenged by the theoretical 
predictions. Efforts are under way to compute the NNLO corrections to top-pair production and it will 
be interesting to compare the ultimately achievable theoretical precision to the experimental accuracy. 
Precise higher order predictions (possibly including resummation), both for inclusive top and top+jets 
production, should also be very valuable for obtaining precise background estimates, such as in Higgs 
searches. Although it will be tried to calibrate the backgrounds with the data themselves, by using 
background-enriched samples for the normalization [43,44], the theoretical predictions are still needed 
for the extrapolation from the background-rich to the signal-enriched regions of phase space. A good 
theoretical precision will lead to reduced systematics on the background, which will be most relevant for 
searches with small signal-to-noise ratios. It is worth mentioning that for the measurement of the b-jet 
cross section similar observations hold as for the top, ie., the statistical error will soon be negligible, 
whereas the systematic uncertainty is expected to be around 15-20%, dominated by the JES. 

Finally, top production will become an extremely valuable calibration tool. The mass peak can 
already be reconstructed with much less than 1 fb _1 , even without b-tagging requirements. With a clean 
sample in hand, it can be exploited for controlling the b-tagging efficiency and serve as a closure test for 
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the JES collections determined from other processes. Concerning the JES, the mass of the hadronically 
decaying W serves as a calibration handle. CMS expects that for intermediate jet pr values this sample 
could lead to JES uncertainties around or below 3% [21]. 

2.5 Conclusions 

I have summarized the experimental and theoretical prospects for some of the most important measure- 
ments of SM processes at the LHC, namely jet, vector boson and top production. The early benchmark 
measurements will include the inclusive jet cross section, the dijet cross section and the dijet ratio, pho- 
ton/Z plus jet production, the Z rapidity distribution, ratios of W and Z cross sections, the Z transverse 
momentum distribution and top pair production. I have indicated the expected uncertainties of the mea- 
surements and shown how these processes serve as tools for the understanding of the detector, for the 
control of backgrounds and for the validation and tuning of Monte Carlo generators. Particularly inter- 
esting are ratio measurements, because otherwise important systematic uncertainties cancel out in this 
case. With differential predictions at higher order in perturbation theory in hand, I have highlighted the 
importance of comparing theory and experiment for the same acceptance cuts, thus avoiding extrapola- 
tion errors. It is important to have (differential) NLO predictions, possibly combined with resummation 
calculations such as implemented in the Monte Carlo generator MC@NLO [45,46], for as many pro- 
cesses as possible. For the cases where this appears to be difficult to achieve, LO plus parton shower 
approaches might still be very valuable tools. However, higher order predictions, up to NNLO, should 
be aimed for as benchmarks, at least in a few cases. I have identified dijet, Z+jet and top production as 
most interesting cases for investing the efforts towards NNLO calculations. 

Acknowledgements 

I would like to thank the organizers for inviting me to this most stimulating workshop, my colleagues 
within the CMS and ATLAS collaborations for many interesting discussions on the topic and Ch. Anas- 
tasiou and M. Dittmar for useful comments on the manuscript. 

Parti 

NEW APPROACHES 

3. ON-SHELL RECURSION RELATIONS 6 
3.1 Introduction 

The efficient calculation of scattering amplitudes with many external legs is a challenging task and 
needed for phenomenological studies at TeV colliders. In the past years, various new methods for 
efficient calculations in QCD have been introduced, originally motivated by the relation of QCD am- 
plitudes to twistor string theory [47]. These methods include the diagrammatic rules of Cachazo, Svrcek 
and Witten (CSW) [48], where tree level QCD amplitudes are constructed from vertices that are off- 
shell continuations of maximal helicity violating (MHV) amplitudes [49], and the recursion relations of 
Britto, Cachazo, Feng and Witten (BCFW) [50,51] that construct scattering amplitudes from on-shell 
amplitudes with external momenta shifted into the complex plane. These developments have triggered 
significant research and numerous applications towards Bom amplitudes in QCD [52-73]. In addition, 
when combined with the unitarity method [74,75] the recursion relations have proven very useful for one- 
loop calculations in QCD [76-105]. Here, we would like to review the basics of the on-shell recursion 
relations for Born QCD amplitudes and the proof of its validity. 

Contributed by: C. Schwinn, S. Weinzierl 
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3.2 Helicity amplitudes and colour decomposition 

It is a well-known fact that the complexity of a calculation based on Feynman diagrams grows factorially 
with the number of external particles. In order to keep the size of intermediate expressions under control, 
a divide-and-conquer strategy has been proven useful: One divides the quantity to be calculated into 
smaller pieces and calculates the small pieces separately. 

One first observes that it is not necessary to square the amplitude and sum over the spins and 
helicities analytically. It is sufficient to do this numerically. This avoids obtaining 0(N 2 ) terms from an 
expression with O(N) terms. The individual amplitudes have to be calculated in a helicity or spin basis. 
This is straightforward for massless fermions. The two-component Weyl spinors provide a convenient 
basis: 

|p±) = i(l±7S)«(p)- (i) 

In the literature there are different notations for Weyl spinors. Apart from the bra-ket-notation there is 
the notation with dotted and un-dotted indices: The relation between the two notations is the following: 

\P+)=PB, (p+\=PAi \P-)=P B , (P~\=P A - (2) 
Spinor products are defined as 

(pq) = ( p - \ q +), [pq] = (p + \q-}, (3) 

and take value in the complex numbers. It was a major break-through, when it was realised that also 
gluon polarisation vectors can be expressed in terms of two-component Weyl spinors [106-112]. The 
polarisation vectors of external gluons can be chosen as 

where k is the momentum of the gluon and q is an arbitrary light-like reference momentum. The depen- 
dence on the arbitrary reference momentum q will drop out in gauge invariant quantities. 

The second observation is related to the fact, that individual helicity amplitudes can be decom- 
posed into group-theoretical factors (carrying the colour structures) multiplied by kinematic functions 
called partial amplitudes [113-117]. These partial amplitudes do not contain any colour information and 
are gauge-invariant objects. In the pure gluonic case tree level amplitudes with n external gluons may be 
written in the form 

A n (l,...,n) = <f~ 2 Yl 2Tr(rMD,..TMn))4( (7 (i) i ... i(7 ( n )) i ( 5 ) 

a&S n /Z n 

where the sum is over all non-cyclic permutations of the external gluon legs and the normalisation of 
the colour matrices is Tr T a T b = 5 ab /2. The quantities A n on the r.h.s. are the partial amplitudes 
and contain the kinematic information. They are colour-ordered, e.g. only diagrams with a particular 
cyclic ordering of the gluons contribute. In general, the colour factors are combinations of open strings 
(T ai ...T an ) i ■_ and closed strings Tr (T bl ...T bm ) of colour matrices. These building blocks form a basis 
in colour space. The choice of the basis for the colour structures is not unique, and several proposals for 
bases can be found in the literature [118-120]. 

3.3 Spinor space versus momentum space 

It will be useful to discuss the relationship between spinor space and complexified momentum space. Let 
us first fix our conventions. The metric tensor is = diag(+l, —1, —1, —1). A null-vector satisfies 

(Pof ~ {Plf ~ (P2) 2 ~ (P3? = 0. (6) 
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This relation holds also for complex p^. In complexified momentum space it is possible to choose a basis 
consisting only of null-vectors: 

ei = (1,0,0,1), e 2 = (0,M,0), e 3 = (0, 1, -i, 0), e 4 = (1, 0, 0, -1), (7) 

is an example of such a basis. Light-cone coordinates are defined as follows: 

P+=P0+P3, P-=P0-P3, P_L=Pl+iP2, P±*=Pl-iP2- (8) 

Note that p±* does not involve a complex conjugation of p\ or p2. A convenient representation for the 
Dirac matrices is the Weyl representation: 



a» 
a* 



), 75 = itVtV = ( J < B = (l,-a), a^ = (l,a), (9) 



with a = (cr x ,a y ,(j z ) being the Pauli matrices. A Weyl spinor pa is an element of a complex two- 
dimensional vector space S, and similar a spinor p B is an element of (another) complex two-dimensional 
vector space S'. We will think of pa and p B as independent quantities. The dual space to S will be 
denoted by S, its elements by p A . Similarly, we denote the dual space to S' by S' and its elements by 
p B . The two-dimensional antisymmetric tensor provides an isomorphism between S and S as well as 
between S' and S': 

p A = e AB PB , pb=P A sab, p A = e AB p B , p B =p A e AB . (10) 
We take the two-dimensional antisymmetric tensor as 



-1 I 



Spinors are solutions of the Dirac equation, therefore we have for massless Weyl spinors 

P^\p+) = 0, p^\p-) = 0, (p+\p^ = 0, (p-\p^ = 0. (12) 
As normalisation we take for massless spinors 

(p ~ Wn\P~) = 2 P^ (P + Mp+) = 2 Pv ( 13 ) 
The solutions to eqs. (12), (13) and (10) are 

e < a -¥) ( - p± » \ _ e -i a +¥j i ^ 



(P+\ = a r i~PX,P+), (p-\ = 7=T (P+>P-L*) • ( 14 ) 

vIp+I vIp+I 

Here a is an arbitrary phase and (f> is the phase of p + = \p + \ e 1 ^. The spinors corresponding to a 
four- vector p^ are only determined up to a phase. With these spinors we have 

(pq)[qp] = 2p-q. (15) 

It is worth to note that the relation u(p) = u(p)^°, or equivalently 

|p+) f = (p+l, \p-) j = (p-\, (16) 
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holds only for real p^ and positive p + (e.g. 4> = 0), since 

|P+) f = „ . *+ *) > \P~) ] = ~ /F=r ( P + * > Pl - ) ■ (17) 

vIp+I vIp+I 

Here the upper asterisk denotes the usual complex conjugation. A pair of spinors {p a ,Pb) determines a 
(unique) null-vector through 

= \pxtf B PB = \(p+\o l ,\p+). (18) 

This is just eq. (13) written reversely. For arbitrary p A and pb the four-vector p^ will be in general 
complex. While eq. (14) defines a map from complexified momentum space to the spinor space S and 
S', which is unique up to a phase, eq. (18) goes in the reverse direction: It defines a map from the space 
S' x S to complexified momentum space. In this context it is worth to observe that if we change p A or 
Pb (but not both) by a linear transformation as 

PA^Pa + Z 1A ° r PB^VB- ZQB, (19) 

the resulting four- vector p^(z) will be a linear function of z. Note however that a linear change in p^ as 
in p^ — > p^ + zq^ with a subsequent application of eq. (14) will not result in a linear change in p A nor 
Pb- 



3.4 On-shell recursion relations 

In the previous section we have seen that we can associate to any null- vector p^ a pair of spinors (p a ,pb)- 
From this pair we can reconstruct the original four-vector through eq. (18). To state the on-shell recursion 
relations it is best not to view the partial amplitude A n as a function of the four-momenta, but to replace 
each four-vector by a pair of two-component Weyl spinors. Therefore the partial amplitude A n , being 
originally a function of the momenta kj and helicities Xj, can equally be viewed as a function of the 
Weyl spinors k° A , k° B and the helicities A.,-: 

A n {k\i Ax, k n , X n ) = A n (k A , k^, Ai, k A , k B , X n ). (20) 

Let us now consider the n-gluon amplitude. For the recursion relation we single out two particles i and 
j. If (A^ Xj) (— , +) we have the following recurrence relation: 

An (k A , k^, Ai, k A , kg, X^nJ = (21) 
[-,k A ,k^,Xi,...,iK A ,iK B ,-Xj j^Ar (^K a ,K b , A, k? A , k J B , Xj, .. 

partitions A=± 

where the sum is over all partitions such that particle i is on the left and particle j is on the right. The 
momentum K is given as the sum over all unshifted momenta of the original external particles, which 
are part of Al. In eq. (21) the shifted spinors k\, k^, Ka and K B are given by 



and 



K 2 
(i+\K\j+) 



(23) 
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Here we shifted k\ and while k 1 ^ and k? B were left untouched. We could equally well have used 

the other choice: Shifting k 1 ^ and k° B , while leaving k\ and k 3 - unmodified. In this case one obtains 
a recursion relation valid for the helicity combinations (A,,, Xj) / (+,—). Therefore for all helicity 
combinations of (Aj, Xj) there is at least one valid recursion relation. Applying this recursion relation to 
the six-gluon amplitude Aq(1~ , 2~, 3~, 4 + , 5 + , 6 + ) with three positive and three negative helicities, we 
choose (i, j) = (6, 1). In this case only two diagrams need to be calculated and we obtain the compact 
result 



4s(l-,2-,3~4 + ,5+ 6+) = 

(6 + |l + 2|3+) 3 (4+|5 + 6|l+) 



4?: 



+ 



[61] [12] (34) (45)5i26 <2 + |1 + 6|5+) [23][34](56)(61)si 56 (2 + |1 + 6|5+) 



(24) 



3.5 Quarks, massive or massless 

QCD does not consist solely of gluons, but contains the quarks as well. Let us now discuss the general 
case of the inclusion of massive quarks. All formulae will have a smooth limit m — ► 0, therefore the 
case of massless quarks will need no further discussion. For massive fermions we have to consider Dirac 
spinors. We can take them as 

"( ± ) = / b-M T~\ ^ + m)\q±), u(±) = -, _! r— (qT\ty + m), 
W =F \q±) (q =F |p>±) 

v(±) = , b \ tt (j>-m)\q±), v(±) = — — ( qT \y-m). (25) 

Here, p is the momentum of the fermion and \q+) and (q + | are two independent Weyl spinors used as 
reference spinors. These two spinors define a light-like four-vector q^ 1 = \{q+ \^ L \q+), which in turn 
is used to associate to any not necessarily light-like four- vector p a light-like four- vector p 9 : 

2 

P b = P-i^Q- (26) 
2p-q 

The reference spinors are related to the quantisation axis of the spin for the fermion, and the individual 
amplitudes with label + or — will therefore refer to this spin axis. From the Dirac spinors we can 
reconstruct the four-vector p^ as follows: 

A 

For the recursion relation, we again single out two particles i and j, which need not be massless, with 
four-momenta pi and pj. To these two four-momenta we associate two light-like four-momenta U and lj 
as follows [121,122]: 



h=, (Pi-oijpj), lj = (-mPi+Pj), a k = -^ — sl g n ( 2 P^j) A _ (28) 

1 - oiotj 1 - OiOtj 2p 2 k 

with A = (2piPj) 2 — 4:pfpj. These light-like four-vectors define massless spinors \k+), (k + |, \lj+) 
and (lj + |. If particle i is a massive quark or anti-quark, we use \lj+) and (lj + | as reference spinors 
for particle i. If particle j is a massive quark or anti-quark, we use \k+) and (Zj + | as reference spinors 
for particle j. We have the recursion relation 

A n (u 1 (-),ui(+),X 1 ,...,u n (-),u n (+),X n ) = (29) 

{■■■,u' i (-),u i (+),X i ,...,iv' K (-),iv K (+),-X) 



partitions A=± 

i 

m 
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Here we denote by k the intermediate particle where we factorise the amplitude, and by K the off-shell 
four-momentum flowing through this propagator in the unshifted amplitude. We shift the Dirac spinors 
as follows: 

K 2 - m 2 

u l '(-)=u i (-)-z\l j +), u' j (+) = u j (+) + z(l t + \, z = jj—^^-j - (30) 

For the intermediate particle k we define the polarisations with respect to the reference spinors \lj+) and 
(h + \: 

= (i fb + \_) W + m k ) \h-) , u' K {+) = {l ] K \, +) + mk ) ' (31) 

where 

K -, = ^_£ (i( + Mi . +)| ^ = ^_l_^_ (i , + | 7 ,| i . +) . (32) 

The recursion relation is valid for (Aj, Xj) / (— , +) with the following exceptions: 

• Particles i and j cannot belong to the same fermion line. 

• The combinations (qf,gf), {qf,gf), (g~,q~) and {gj ,qj) are excluded. 

• If i is massive, the combinations {qf , q'j + ), {qf, <7j + )> {qf , q'j + ) and {qf,q'j + ) are excluded. 

• If j is massive, the combinations (q~ , q'j~), {q^,q'j~), (q^,q'j~) and (q~,q'j~) are excluded. 
Instead of shifting Ui{— ) and Uj(+), we can alternatively shift Ui(+) and Uj(-): 

K 2 - m 2 

u' i {+) = u i {+)-z(l j + \, u/{-)=u j {-) + z\l i +), z= y (33) 

For the intermediate particle k we define the polarisations now with respect to the reference spinors \k+) 
and (lj + | : 

uk\-)= {^b + ij ) {H' + m k )\l ] -), u' K {+) = _ 1 (U-\ (#' + m k ) , (34) 

where 

K" = K»-l(l i + m i i+) , ^ = ^-\ {h+ f m+) (h + W>\k+). (35) 

Doing so, we obtain a recursion relation valid for (A;, Xj) ^ (+, — ) with the following exceptions: 

• Particles i and j cannot belong to the same fermion line. 

• The combinations (gf,qf), {gf,qf), [q[-,g]) and {q~,gj) are excluded. 

• If j is massive, the combinations {qf , q'j + ), {qf , q'j + )> {qf , q'j + ) and (qf ,q'j + ) are excluded. 

• If % is massive, the combinations (q~ , q'j~), {q^ , <7j~)> {<li , 1j~) an( ^ {oT i<fj~) 316 excluded. 

As we are free to choose the particles i and j, we can compute all Born helicity amplitudes in QCD 
with two-particle shifts via recursion relations, except the ones which involve only massive quarks or 
anti-quarks. Amplitudes consisting solely of massive quarks and anti-quarks and with more than six 
particles may be calculated recursively if one allows more general shifts, where more than two particles 
are shifted. 
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3.6 Proof of the on-shell recursion relations 

For the proof [51,55,59-61,72] of the on-shell recursion relation we discuss as an example the case of 
the holomorphic shift as in eq. (22) or eq. (30). One considers the function 

A(z) = A n (...,Ui(-),Ui(+),X i ,...,Uj(-),Uj(+),Xj,...) (36) 
of one variable z, where the z-dependence enters through 

Ui(-) = Ui{-) -z\lj+), u'j(+) = Uj(+) + z(U + \. (37) 

The function A(z) is a rational function of z, which has only simple poles in z. This follows from the 
Feynman rules and the factorisation properties of amplitudes. Therefore, if A(z) vanishes for z — > oo, 
A{z) is given by Cauchy's theorem as the sum over its residues. This is just the right hand side of 
the recursion relation. The essential ingredient for the proof is the vanishing of A(z) at z — ► oo. If 
(Aj, Xj) = (+, — ) it can be shown that each individual Feynman diagram vanishes for z — > oo. For the 
helicity combinations (+, +) and (— , — ) one first constructs a supplementary recursion relation based on 
three-particle shifts and deduces from this representation the large z-behaviour of A(z). This establishes 
the recursion relation for these helicity combinations with the exceptions indicated above. The proof for 
the anti-holomorphic shift as in eq. (33) proceeds analogously. 

4. ON-SHELL RECURSION TO DETERMINE RATIONAL TERMS 7 

On-shell methods offer an auspicious approach for dealing with the rapid growth in complexity of loop 
amplitudes as the number of particles in the process increases. These methods rely on the unitarity of 
the theory [123, 124] which requires that the poles and branch cuts of amplitudes correspond to physical 
propagation of particles. On-shell methods are presently undergoing intense development for use at loop 
level (see, for example, refs. [86,87,89,95,97-104, 125-127]). Their advantage lies in the relatively mild 
growth in complexity as the number of external particles increases, effectively reducing loop calculations 
to tree-like calculations. 

On-shell methods fall into two basic categories: the unitarity method [74, 128] which constructs 
amplitudes based on their branch cuts, and on-shell recursion [50, 5 1 ] which constructs amplitudes from 
their poles. In this section we discuss using on-shell recursion as a means for computing rational terms 
of one-loop amplitudes. The loop-level construction is based directly on the construction of tree-level 
recursion relations by Britto, Cachazo, Feng and Witten (BCFW), though a number of new features are 
present. Further discussion of the unitarity method approach, as well as other new methods exploiting 
on-shell conditions on intermediate states [80, 89,95,98, 102, 104, 129] may be found in other sections 
of this report. Introductions to on-shell methods may be found in various reviews [130-132]. Earlier 
reviews of spinor methods, which are profitably used in conjunction with on-shell methods, may be 
found in refs. [133,134]. 

In the context of the unitarity method, it is convenient to divide the amplitudes into pieces that 
contain branch cuts, plus rational (non-cut-containing) pieces. When using dimensional regularization, 
the branch-cut containing pieces may be computed by ignoring the distinction between D = 4 — 2e 
dimensions and four-dimensions in the numerators of the loop-momentum integrands [74, 128]. This 
observation allows powerful four-dimensional spinor techniques to be used to greatly simplify the on- 
shell tree amplitudes appearing in the unitarity cuts. However, if one wants to obtain also the rational 
terms directly from the cuts [130, 135], then the (— 2e) dimensional contributions are crucial: dropping 
these pieces leaves undetermined additive rational terms. (The branch cuts can determine rational terms 
at 0(e) because they develop branch cuts at 0(e).) By using amplitudes valid in D = 4 — 2e dimensions 
in the unitarity cuts, all rational terms are kept 8 but at the cost of more complicated expressions. It has 

'Contributed by: Z. Bern, L.J. Dixon 

8 In the language of dispersion relations [136, 137], this reconstruction is possible because the dispersion integrals converge 
with dimensional regularization [138]. 
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physical poles 



1*1 



spurious poles 




Fig. 6: Using Cauchy's theorem, rational terms in loop amplitudes can be reconstructed from residues at poles in the complex 
plane. The poles are of two types: physical and spurious. All pole locations are known a priori. Residues at physical poles 
follow from factorization onto lower-point amplitudes. Residues at spurious poles cancel against corresponding contributions 
from the cut parts, and so they can be inferred from four-dimensional cuts. 

been pointed out [90-93] that the rational terms are relatively easy to obtain from Feynman diagrams 
because they do not require the full set of tensor integrals. In addition, Brandhuber et. al. have argued 
that the rational terms can be obtained from a set of counterterms [125]. Britto and Feng have recently 
given a complete set of formulae for constructing loop amplitudes, including their rational terms [104], 
following earlier work [80, 89, 97-100, 129]. 

An early version of on-shell methods was used to compute the one-loop matrix elements needed 
for the NLO QCD corrections to e+e~ — ► 7*, Z — > 4 jets and pp — ► W, Z + 2 jets [139]. They have also 
been used to obtain analytic expressions for the complete one-loop six-gluon amplitude [74,76,78, 80, 
86,87,89-92, 128] as well as a variety of helicity configurations for n-gluon amplitudes [82,85-87, 101]. 
The results confirm the mild growth in complexity of these methods as the number of external particles 
grows. 

A crucial next step for applying these methods to LHC physics is the construction of automated 
programs to compute the large number of phenomenologically interesting high-multiplicity processes. 
As discussed in other sections of this report, such automated programs are in the midst of being con- 
structed [103, 126], using the integration machinery of Ossola, Papadopoulos, and Pittau [95]. The recent 
numerical implementation by Ellis, Giele and Kunszt [126] of the unitarity method presently makes use 
of D = 4 simplifications and hence does not contain rational terms. The program of Ossola, Papadopou- 
los, and Pittau [103] can be used to obtain the rational terms, but currently requires one-loop Feynman 
diagrams to capture these terms, instead of more efficient on-shell tree amplitudes. 

On-shell recursion offers an efficient alternative for constructing one-loop rational terms directly 
from their known factorization properties, in much the same way as the BCFW recursion relations can 
be used to obtain tree-level amplitudes. However, a number of new issues arise at loop level that must 
be dealt with first to have a practical method. These issues include the appearance of branch cuts, spuri- 
ous singularities, the behavior of loop amplitudes under large complex deformations and in some cases, 
'unreal poles', which are present with complex but not real momenta. More practical issues are automa- 
tion and numerical stability. Here we briefly summarize the construction of rational terms via on-shell 
recursion [81-83, 86, 87, 101], describing in particular a simple modification making it straightforward 
to automate. 

In general, any one-loop amplitude can be divided into two pieces, 



where C n are the 'cut-containing terms' possessing logarithms, polylogarithms, and associated it's. The 



n 



c r C n + R 



-11 



(38) 
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Fig. 7: Diagrammatic contributions to on-shell recursion at one-loop for a [j, I) shift. The labels 'T and 'L' refer to tree and 
one-loop vertices corresponding to the rational parts of lower-point on-shell amplitudes. 



rational terms, denoted by R n , are denned by setting these (poly)logarithmic terms to zero, 



cr 



= A n 
rat Cr 



(39) 



ln,Li2 ,7r— >0 



Let us assume that the cut-containing terms C n of the particular amplitude under consideration have 
already been computed using four-dimensional unitarity. This leaves the problem of computing the 
rational terms R n . 

On-shell recursion relations can be derived by considering complex on-shell deformations of am- 
plitudes A{z), which are characterized by a single complex parameter z [51]. The z-dependence al- 
lows us to use standard complex variable theory to construct amplitudes via Cauchy's Theorem. To set 
up an on-shell recursion relation for R n consider the effect of shifting some set of external momenta 
h — > ki(z), such that the on-shell conditions [/cj(z)] 2 = mf and the original momentum conservation 
are satisfied. In the massless case, it is particularly convenient to shift the momenta of two external legs, 
say, j and I, 



;{r\r\n, 



K{z) = ^ + Uj-\^\n, 



(40) 



where z is a complex parameter and \i + ) and \i ) are Weyl spinors of positive and negative chirality, 
following the notation of ref. [133]. In terms of these spinors, the shift is 



i r >-i r >-*|r> 



(41) 



We denote the shift in eqs. (40) and (41) as a [j, I) shift. 

The on-shell recursion relations follow from evaluating the contour integral, 



1 

2tH 



dz 



Rn(z) 



(42) 



c 



where the contour is taken around the circle at infinity, as depicted in fig. 6, and R n (z) is R n evaluated 
at the shifted momenta (40). If the rational terms under consideration vanish as z — > oo, the contour 
integral vanishes and we obtain a relationship between the desired rational contributions at z = 0, and a 
sum over residues of the poles of R n (z), located at z a , 



Rn{0) = ~ ReS * 
poles a 



Rn{z) 



(43) 



If R n {z) does not vanish as z — > oo, then there are additional contributions. A systematic strategy for 
computing such large z contributions using auxiliary recursion relations was presented in ref. [86], to 
which we refer the reader. 
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Fig. 8: Diagrammatic representation of one-loop overlap terms for a [j, I) shift. The channels correspond to physical poles and 
remove double counts induced by cut completion. 



As illustrated in fig. 6, the poles in one-loop rational terms fall into two categories: the physical 
poles, which are present in the full amplitudes; and the spurious poles, which cancel against poles in the 
cut-containing terms. 

Residues at physical poles are dictated by factorization onto lower-point amplitudes. They may be 
computed using the recursive diagrams 9 in fig. 7, 



phys. poles {r,s} 



£ I A%»(z = z rs ) -4- Rr(z = z rs ) + R L {z = z rs ) -L_ A^{z = z 



1Y r---s r---s 



+ AT c (z = z rs ) ^f- AT(z = z rs ) \ . (44) 
r—s \ 

The 'vertices' Ri and Rr in this recursion relation are the pure rational parts — using the definition (39) 
— of the lower-point, on-shell one-loop amplitudes. The 'vertices' A^ 00 and ^4^ ee are on-shell tree 
amplitudes. The subscripts L and R on the vertices indicate their location to the left or right of the 
central propagator in fig. 7. In the vertices the shift variable z is frozen to the values 

K 2 

r " s (45) 



(rw r ... s \i-)' 

corresponding to the location of the poles in z, coming from shifted propagators. The rational part Rjr 
of the factorization function T [140] only contributes in multi-particle channels, and only if the tree 
amplitude contains a pole in that channel. Generically we have a double sum, labeled by r, s, over 
recursive diagrams, with legs j and I always appearing on opposite sides of the pole. There is also a sum 
over the helicity h of the intermediate state. The superscript D on R® indicates that this set of recursive 
diagrammatic contributions is not the whole rational part, as discussed below. 

It is interesting to note the similarity of the one-loop recursion relation (44), to the corresponding 
tree-level recursion relation [51], 

r,s,h r - s 

Thus loop-level recursive diagrams echo the simplicity of tree-level recursion. 

One way to deal with the spurious poles is to start by finding a 'cut completion' C n [83, 86, 87, 
101, 132]. One adds certain rational terms CR n to C n , such that the spurious poles in C n {z) cancel 
entirely. Because physical amplitudes cannot have spurious singularities, the remaining rational terms, 



''Unreal' poles, which do not correspond to factorizations with real momenta, may be avoided by choosing appropriate 
shifts [86]. 
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Fig. 9: Non-vanishing recursive diagrams for the rational terms of A^' a (l , 2 , 3 + , 4 + , 5 + ), using a [1, 2} shift. 



Rn(z), must also be free of these spurious singularities. This cut completion makes it unnecessary to 
compute residues at spurious poles. It is rather helpful when deriving compact analytic expressions for 
the amplitudes. It does introduce additional 'overlap diagrams', as depicted in fig. 8. These diagrams 
correct for the contributions of CR n in physical factorization limits. They are simple to compute from 
the residue of CR n at each physical pole z rs . 

Following the cut-completion procedure, a variety of rational terms with an arbitrary number of 
external legs have been constructed [85-87, 101], giving complete amplitudes when combined with the 
previously-computed cut-containing parts [74,78,79,84, 101, 128]. More generally, it should be possible 
to form a set of cut completions using integral functions of the type given in ref. [141] to absorb the 
spurious singularities. 

For the purposes of automation in a numerical program, another approach is preferable [142]. It 
is simpler to obtain the residues at the spurious poles directly from the cut parts, calculated from the 
four-dimensional unitarity method. Because a complete amplitude is free of spurious poles, any spurious 
pole found in the rational parts must cancel a spurious pole in the cut parts. To get the full rational part, 

Rn = R n + R n , (47) 

we add to the recursive diagrams R® some 'spurious' contributions Rf t , evaluated by means of the cut 
terms C n (z), 

RS = - £ Res_^= £ Res z=z ^. (48) 

spur, poles /3 spur, poles /3 

The spurious poles j3 can be classified systematically in terms of the vanishing loci, A(z) = 0, of shifted 
Gram determinants A associated with box, triangle and bubble functions. (In the massless case, the 
bubble Gram determinant does not generate a spurious pole.) 

To illustrate this modified procedure, consider the five-gluon amplitude A^' s (l~ , 2~, 3 + , 4 + , 5 + ), 
with a scalar in the loop. The construction of the rational terms in this amplitude, using on-shell recursion 
with cut completion, has already been discussed in some detail elsewhere [83, 132]. Here we describe 
the new approach for obtaining these terms. 

The cut part of the amplitude [143] is 
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where /i 2 is a scale and ' • • • ' signifies that we are dropping terms not pertinent for our discussion. The 
spinor inner products and kinematic invariants are defined as, 

(ab) = (a~ \b + ) , [ab] = (a+\b-), s ab = (k a + k b ) 2 . (50) 
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The rational terms are determined by evaluating the recursive diagrams, plus the rational residues 
of the cut terms at the spurious poles. Here we use the [1, 2) shift. (As discussed in ref. [83, 86], for 
this shift there are no additional contributions from either large z behavior or unreal poles.) With the 
[1, 2) shift, the non- vanishing recursive diagrams are depicted in fig. 9. A simple computation of these 
diagrams (see section 5.1 of ref. [132]) gives, 



5 \3e 9/ (23)(34)(45)(51) ' 5 3 (3 4)[1 2][1 5] [2 3] 2 ' V ' 

as the recursive contributions. 

We still need to account for the residues at the spurious poles. In the present example with a [1,2} 
shift, the only such pole comes from solving s^i(z) — S2s(z) = (corresponding to a shifted two-mass 
triangle Gram determinant). The solution is, 

_ S 5 l - S 2 3 S 5 l - S 2 3 , 52 - 



A (15)[5 2] + (13>[3 2] (14)[24] 

To obtain the residue, we start from the logarithmic terms of eq. (49), and perform the [1, 2) shift eq. (41), 
yielding, 

C&(z) _ i [34](41)((24)+z(14))[45](((23)+z(13))[34](41)+((24)+z(14))[45](51)) 
z 3 (34) (4 5) 

, «2 3)+z(13))[3 2] 
1 <5l)([l5j- Z [25j) 

+ • • • , (53) 
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where we have kept only the term contributing to the spurious residue at z s . 

The residue needed for eq. (48) can be extracted straightforwardly, by series expanding both the 
logarithm and its coefficient in eq. (53) around z = z s . Cleaning up the result of this residue evaluation, 
we find, 

4 a) = Res^^M 
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The total rational part, 



(54) 



R 5 = Rl? + Rl = + d£ ] + , (55) 

matches the result obtained in refs. [83, 132] using a cut completion. The complete amplitude is obtained 
by summing the cut (49) and rational (55) contributions. 

The modified construction described here is amenable to automation. In a numerical program, 
instead of obtaining the residues at spurious poles by series expansion, we may compute them by nu- 
merically evaluating the cut terms at several points around each pole. The automation and numerical im- 
plementation of on-shell recursion to amplitudes of interest for LHC phenomenology will be described 
elsewhere [142]. 



Acknowledgments 

We are grateful to Carola Berger, Fernando Febres Cordero, Darren Forde, Harald Ita, David Kosower 
and Daniel Maitre for collaboration on the topics discussed here. 



28 



5. FOUR- AND D-DIMENSIONAL UNITARITY CUTS 



5.1 Four-dimensional unitarity cuts 

The application of unitarity as an on-shell method of calculation, as introduced in [74], is based on the 
principles that products of on-shell tree-level amplitudes produce functions with the correct branch cuts 
in all channels [144-147], and that any one-loop amplitude can be expressed as a linear combination of of 
scalar (i.e. trivial numerator) master integrals [148-153]. Given the independent knowledge of the master 
integrals, to compute any amplitude it is sufficient to evaluate the coefficients of such a decomposition. 

For one-loop amplitudes, systematic techniques have been developed to extract the coefficients 
algebraically, preserving gauge invariance at every intermediate stage of the computation. The use of 
four-dimensional states and momenta in the cuts enables the construction of the polylogarithmic terms 
in the amplitudes, which are fixed by their branch cuts, but generically drops rational terms, which have 
to be recovered independently. 

Some recent developments of unitarity-based methods apply generalized unitarity cuts to ampli- 
tudes and master integrals. The coefficients are then extracted by matching the generalized cuts. General- 
ized unitarity corresponds to requiring more than two internal particles to be on-shell, and the fulfillment 
of these constraints can only be realized through complex kinematics. Complex kinematics are the key 
for the exploration of singularities of amplitudes and the use of factorization information to reconstruct 
amplitudes recursively, since the singularities of a scattering amplitude are determined by lower-point 
amplitudes in the case of poles and by lower-loop ones in the case of cuts [49, 130, 133, 134]. 

A notable application of complex momenta within generalized unitarity is the quadruple cut, which 
allows for an immediate and purely algebraic determination of the coefficients of box functions [129]. 
Every box coefficient is simply determined by the product of the four tree-level amplitudes sitting at 
each corner, evaluated at the two particular values of the loop momentum which fulfill the four equa- 
tions imposed by the vanishing of the cut denominators. Double and triple unitarity cuts have led to 
direct techniques for extracting triangle and bubble integral coefficients analytically [80,89,99]. In cases 
where fewer than four denominators are cut, the loop momentum is not frozen, so some explicit integra- 
tion over the phase space is still required. In [80, 89,99], double or triple cut phase-space integration 
has been reduced to extraction of residues in spinor variables, and, in the case of a triple cut, residues 
in a Feynman parameter. This approach has been used to compute analytically the final contributions 
to the cut-constructible part of the the six-gluon amplitude [80, 89], and the complete six-photon ampli- 
tudes [93,94]. 

In general, one can compute n-point (n > 4) coefficients from quadruple cuts, three-point coeffi- 
cients from triple-cuts, and two-point coefficients from double-cuts, by avoiding the conventional tensor 
reduction. As it turns out, given the decomposition of any amplitude in terms of master integrals, the 
coefficient of any n-point master integral can be recovered from the ?7,-particle cut. Obviously, any n- 
particle cut may also detect higher-point master integrals, which appear with different analytic structures 
for they come from the Landau poles specific to each of the master integrals. This is indeed the case 
for the usual (double) unitarity cut, which can be used exclusively to derive box, triangle, and bubble 
coefficients. In cases with massive particles, it is useful to apply a generalized cut to find the coefficient 
of the 1 -point (tadpole) master integral. 

The algorithm of [80, 89] for evaluating any finite unitarity cut involves a change of coordinates 
that brings the loop momentum variable into the spinor formalism. The idea is that the final integrals 
always localize to some poles in the region of integration. Phase space integration is thus reduced to 
a sequence of algebraic manipulations, up to an integration over a single Feynman parameter, which is 
responsible for logarithms. Ultimately, even this integration does not need to be carried out, since it 
is possible to match integrands at an early stage of the calculation. The procedure naturally leads to 
a clean separation of the master integrals, allowing for an individual calculation of the corresponding 
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coefficients. 

By now, explicit analytic formulas for the results of unitarity-based methods are available [100, 
102, 104, 105, 129]. Coefficients of the master integrals are listed directly in terms of tree-level input 
data. All integration and reduction can now be avoided. Although it may not be a significant distinction 
in terms of the final results, we note that the derivations of [102, 105, 129] used generalized cuts, while 
those of [100, 104] used ordinary double cuts. 

5.2 D-dimensional unitarity 

Full one-loop amplitudes can be reconstructed from unitarity cuts in D = 4 — 2e dimensions [135, 138]. 
In the D-dimensional unitarity method, there is no need to distinguish "rational" and "cut-constructible" 
parts of the amplitude. Contributions that might be called "rational" (after expanding around e = 0) 
appear here as e-dependent terms in the coefficients of the master integrals (before expanding around 
e = 0). 

A systematic D-dimensional unitarity double-cut method was proposed in [97,98], reducing one- 
loop amplitudes to master integrals for arbitrary values of the dimension parameter. Coefficients of the 
master integrals can be extracted without fully carrying out the D-dimensional phase space integrals. 
Only a four dimensional (massive) integration is explicitly required. That can be performed by four- 
dimensional unitarity techniques or any other available alternative. The remaining integral, which gives 
rise to the e-dependence of the cut-amplitude, is mapped to phase-space integrals in 4 + 2n — 2e dimen- 
sions, where n is a positive integer. With recursive dimensional shift identities, similar to the ones in 
loop integration, the cut-amplitude is reduced in terms of bubble, triangle, box and pentagon cut master 
integrals in 4 — 2e dimensions. The reduction is valid for an arbitrary number of dimensions. Expand- 
ing in e gives both the (poly)logarithmic and rational part of the amplitude at O(e ) and higher; these 
contributions are required in cross-sections beyond the next-to-leading order in the relevant coupling 
strength. 

Generalized unitarity cuts are possible and useful in D dimensions as well [99, 154]. The benefits 
of the double-cut integration of [80, 89,97,98] have been extended to the evaluation of triple cuts [99], 
for the direct extraction of triangle and higher-point-function coefficients from any one-loop amplitude 
in arbitrary dimensions. Accordingly, the triple cut is treated as a difference of two double cuts with the 
same particle content, and the same propagator carrying respectively causal and anti-causal prescription 
in each of the two cuts. The triple cut phase space for a massless particle in D dimensions is written 
as a convolution of a four-dimensional triple cut of a massive particle, and an integration over the cor- 
responding mass parameter, which plays the role of a (— 2e)-dimensional scale. Just as in the case of 
the double-cut [97,98], to perform the four-dimensional integration, one combines the method of spinor 
integration of massive phase-space integrals, and an integration over the Feynman parameter. But, in the 
case of the triple-cut, after Feynman parametrization, by combining back the two double-cuts, the para- 
metric integration is reduced to the extraction of residues to the branch points in correspondence of the 
zeroes of a standard quadratic function in the Feynman parameter. It is that standard quadratic function 
(or rather, its roots) that carry the analytic information characterizing each master integral, therefore de- 
termining its own generalized cuts. The final integration over the dimensional scale parameter is mapped 
directly to the triple cut of master integrals, possibly with shifted dimensions. 

5.3 Mathematica package for spinor formalism 

Recently, the package S@M (Spinors® Mathematica) was released [155]. It implements the spinor- 
helicity formalism in Mathematica. The package allows the use of complex-spinor algebra along with the 
multi-purpose features of Mathematica, and it is suitable for the algebraic manipulation and integration 
of products of tree amplitudes with complex spinors sewn in generalized unitarity cuts. 
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6. COMMENTS ON UNITARIT Y BASED ONE-LOOP ALGORITHMS 1 1 



6.1 Introduction 

At the LHC deviations from the Standard Model will likely show up in observables of complex multi- 
particle final states. It is important to understand the Standard Model predictions and uncertainties for 
these complicated final states. Leading-order Monte Carlo (LO-MC) programs give a first estimate. 
However, to understand the uncertainties we need at least a next-to-leading order Monte Carlo (NLO- 
MC). 

The basic calculational framework for both tree-level amplitudes (needed for the LO-MC) and 
one-loop amplitudes (needed for NLO-MC) is the perturbative expansion in Feynman diagrams. This 
immediately gives us a straightforward algorithm suitable for numerical implementation. However, such 
implementations are not satisfactory from a numerical standpoint. The number of Feynman diagrams 
grows faster than factorial with the number of external particles involved in the scattering process. As a 
consequence the number of multiplications, and therefore the computer time needed to evaluate a phase 
space point, will grow at least as fast. 

In computer science, algorithms with factorial growth are called exponential or factorial algorithms 
or simply E-algorithms [156]. Such algorithms are not considered optimal, i.e. the number of external 
particles we can calculate becomes quickly limited by computer resources. In contrast, the other class 
of algorithms with polynomial growth in the number of external legs are called P-algorithms. Such 
algorithms are highly desirable as the added computational effort needed to go from N to (iV+1) external 
particles is (^r") ■ This means the limiting factor for these types of algorithms in scattering amplitude 
calculations is often human resources instead of computer resources. In the subsequent sections we will 
argue that for numerical solutions, especially in the era of LHC physics, the complexity of the algorithms 
are an important consideration. 

6.2 Tree-level algorithms of polynomial complexity 

The number of Feynman graphs grows very fast with the number of external legs. For a tree-level 
A-gluon scattering the number of individual Feynman graphs is approximately N^ N ~ 3 ^ (within 5% 
accuracy up to 16 gluons) [157]. This means that to extend the LO-MC from 2 gluon — > 5 gluon to 
2 gluon — > 6 gluon, the number of multiplications increases by at least a factor of 13. Several LO-MC 
are available for the numerical evaluation of arbitrary tree-level processes in the Standard Model and 
some of its extensions. Most of these packages are based on simple Feynman diagram evaluations. We 
call these Numerically Implemented Exponential (NIE) algorithms. A prominent representative in this 
class of algorithms is MadGraph [158]. 

By using currents instead of amplitudes in Feynman diagram calculations one can construct re- 
cursion relations [159]. This method re-uses recurring groups of off-shell Feynman graphs in an optimal 
manner. Because this leads to a more factorized way of calculating the scattering amplitude one can im- 
mediately extend the analytic calculations to more complex processes such as vector boson production 
with up to 6 partons [160, 161] and 7 partem processes [162]. 

Another consequence of the recursion relations is the formulation of an algorithm of polynomial 
complexity. For a tree-level A^-gluon process the number of multiplications grows as N A [157]. This 
means that to extend the LO-MC from 2 gluon — > 5 gluon to 2 gluon — > 6 gluon the increase in 
the number of multiplications is only 1.7 (compared to 13 for standard Feynman graph calculations). 
We will denote the LO-MC programs based on recursive type of evaluation Numerically Implemented 
Polynomial (or NIP) algorithms. A prominent representative is the ALPGEN program [163]. 

As is clear from the discussion we have reached a point for LO-MC where the problem of numer- 
ically calculating the scattering amplitudes can be considered solved. 
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6.3 Toward one-loop algorithms of polynomial complexity 

The LO-MC prediction at LHC type of energies for QCD and/or Electro-Weak processes are rather 
qualitative. One estimates the magnitude of the cross section and predicts the shape for an observables. 
The NLO-MC will give us a first real estimate of the expected normalization and will give an order as 
correction to the shape. Within the perturbative context this allows us to estimate the uncertainties on the 
predictions with some confidence. 

The one-loop amplitude of the basic 2 gluon — > 2 gluon was already calculated analytically in 
1986 [164] using the standard Feynman diagram calculation. One can extend this method brute force 
with modern day computers. Using a combination of e.g. QGRAF [165] and FORM [166] one can 
generate and manipulate the Feynman graphs giving tensor coefficients times tensor integrals. The tensor 
integrals can be determined using Pasasarino-Veltman reduction [167] or other techniques. This then 
can be straightforwardly implemented in a numerical code for e.g. 2 gluon — ► 4 gluon [168]. The 
evaluation of a single phase space point for this process is of the order 9 second (10,000 times slower as 
the 2 gluon — ► 2 gluon one-loop amplitude generated using the same procedure). It is clear that such a 
direct approach using Feynman diagrams is severely affected by the factorial growth in complexity. One 
needs badly a polynomial complexity calculational approach. 

It can be shown that any dimensional regulated multi-loop amplitude is fully reconstructible using 
unitarity cuts [169]. Because the unitarity cuts factorizes the one-loop amplitudes into a product of 
two tree-level amplitudes this proves the existence of a polynomial complexity algorithm for one-loop 
calculations. This was exploited in the analytic calculation of the e + e~ — ► 4 partons one-loop amplitude 
[139] 12 . The method applies four-dimensional unitarity cuts, thereby it only partly reconstructs the one- 
loop amplitude through unitarity, the so-called cut-constructible part. The missing part is referred to as 
the rational part and is determined by other methods. The applied 4-dimensional unitarity method has no 
direct numerical equivalent, but it is explicitly demonstrated that such methods of polynomial complexity 
work very well within the context of analytic multi-leg one-loop calculations. 

The first numerically implementable method came from the so-called quadruple cut method [129]. 
While presented as an analytic method to calculate coefficients of the 4-point scalar master integrals for 
multi-gluon processes, it has a direct numerical implementation. The numerical procedure can be used 
to calculate the box coefficients for any multi-particle scattering process. By applying the quadruple 
cut the one-loop graph breaks down into four tree-level amplitudes. This is therefore instantly a NIP 
algorithm for calculating the coefficients of the 4-point master integrals. From the unitarity constraints, 
i.e. the four cut propagators have to be numerically zero, one gets only two complex solutions for the 
loop-momentum. By evaluating the product of the four tree-level graphs using the two complex loop 
momenta solutions, one gets the coefficient by simply averaging over the two solutions. 

The numerical implementation of the method is extremely fast and simple, showing the potential of 
a full numerical implementable unitarity method. To achieve this one also has to calculate the coefficients 
of the other 3 master integrals (the 1-, 2- and 3-point scalar integrals). A direct generalization of the 
quadruple cut method becomes complicated because of overlapping contributions. By applying a triple 
cut to determine the 3-point coefficient one has to take into account that part of this contribution is also in 
the quadruple cut. Disentangling these overlapping contributions proves to be not that straightforward. 

For a one-loop amplitude one can construct a general parametric form of the integrand and deter- 
mine its coefficients by demanding different combinations of sets of propagators to be zero (i.e. cutting 
the lines) for both the parametric form of the amplitude and the expression obtained using Feynman 
graphs [95]. This method is purely algebraic as it works on the integrand level. When setting four propa- 
gators to zero this method is identical to the quadruple cut method. However, we now get in addition the 
full loop dependence of the integrand of the 4-point master functions through its parametric form. This 
allows one to simply determine the triple cut contribution of the parametric 4-point integrand and hence 

l2 The 5 gluon one-loop was calculated using string inspired methods [143]. 
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we know the subtraction term. 

Using this method to construct the subtraction terms it is now straightforward to formulate a nu- 
merical implementable algorithm of polynomial complexity for the cut constructible part [126]. Because 
we determine the coefficients of the 2-, 3- and 4-point parametric form of integrands by the equivalent of 
unitarity cuts, the actual one-loop amplitude factorizes in a product of two, three or four tree-level ampli- 
tudes. That is, we can determine the full parametric form of the integrand from tree-level amplitudes. The 
final loop integration over the parametric form is straightforward and gives us the three scalar master in- 
tegrals and their respective coefficients. This method now extends the polynomial complexity algorithm 
of the quadruple cut method to include also the triple and double cut contributions. As a demonstra- 
tion we used this method to numerically evaluate multi-gluon scattering amplitudes. We found using a 
single standard processor the following results: the 2 gluon — ► 2 gluon at 9 seconds/10,000 events, the 
2 gluon — ► 3 gluon at 35 seconds/10,000 events and the 2 gluon — > 4 gluon at 107 seconds/10,000 events. 
This can be approximated by JV 6 /450 seconds/10,000 events, which by extrapolation would give around 
260 seconds/10,000 events for 2 gluon — ► 5 gluon. These evaluation times are more than sufficient for 
use in NLO-MC generators, even on a modest single processor system. 

6.4 Conclusions: the rational part 

The final step is a numerical suitable algorithm for the rational part of the one-loop amplitude. This is 
the final hurdle in achieving a full solution of polynomial complexity for numerical one-loop amplitude 
evaluations. Three methods exist in the literature. The first method determines the rational part of the 
tensor integrals. These rational parts can then be contracted in with the tensor coefficients to give the full 
one-loop rational part [90,93]. This method goes back to the Feynman diagram expansion and leads to 
an algorithm of factorial complexity. This negates all progress made with the determination of the cut 
constructible part using numerical unitarity techniques. 

The other two methods are more analytic in concept, but should in principle be suitable for a nu- 
merical implementation. The so-called bootstrap method sets up a recursive procedure for the rational 
part [86] similar to the tree-level unitarity based recursion relations [51]. This makes the method of poly- 
nomial complexity. However, in its current formulation it is not suitable for numerical implementation. 
The reason is that both the rational and cut constructible part of the one-loop amplitude contain so-called 
spurious poles. When adding the two parts together these spurious poles cancel. This means that for 
the construction of an unitarity based recursion relation in the rational part these spurious poles have to 
be removed. This procedure is called cut-completion, i.e. make both cut-constructible and rational part 
free of spurious poles. Then the rational part contains only physical poles and a unitarity based tree-level 
like recurrence relation for the rational part is constructible. Unfortunately the cut-completion procedure 
requires analytic knowledge of the spurious terms, which up to now have only be determined by explicit 
analytic calculation of the cut-constructible part. 

One can in principle retrieve the full one-loop amplitude by applying D-dimensional unitarity 
cuts [98, 135]. Such an implementation is per construction of polynomial complexity. It requires the 
calculation of the D-dimensional tree-level amplitudes. This can be implemented by restricting oneself 
to massive scalar internal particles where the mass in generated by the extra-dimensional length of the 
loop-momentum. In this manner the extra-dimensional part of the loop-momentum can be integrated out. 
After that one can read off the appropriate master integral coefficients and rational part. The required 
scalar internal particles restrict this method at the moment to purely gluonic scattering amplitudes. In 
its current implementation this method is restricted to analytic applications for purely gluonic one-loop 
scattering amplitudes. 

It is clear from the discussions that a numerical algorithm of polynomial complexity is the only 
issue left in fully solving one-loop calculations in a similar way tree-level calculations have been solved. 
Achieving this final step would open the way to a multitude of NLO-MC generators for processes such 
as for example PP — ► tt +2 jets, PP — ► tt + bb and PP — ► Vector-Boson + 3,4 jets. 
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A NIP implementation for the rational part has to exist. Its construction in the near future is of 
great importance to make the first step towards more complicated NLO-MC programs relevant for the 
LHC phenomenology. 



7. PHYSICAL APPLICATIONS OF THE OPP METHOD TO COMPUTE ONE-LOOP AM- 
PLITUDES 13 

7.1 Introduction 

In two recent papers [95, 96], we proposed a reduction technique (OPP) for arbitrary one-loop sub- 
amplitudes at the integrand level [121] by exploiting numerically the set of kinematical equations for 
the integration momentum, that extend the quadruple, triple and double cuts used in the unitarity-cut 
method [105, 126, 128, 129]. The method requires a minimal information about the form of the one-loop 
(sub-)amplitude and therefore it is well suited for a numerical implementation. The method works for 
any set of internal and/or external masses, so that one is able to study the full electroweak model, without 
being limited to massless theories. 

In Section 7.2 we outline the basics features of the method. In Section 7.3 we describe a numeri- 
cally stable implementation of the OPP algorithm, in a form of a FORTRAN 90 code, CutTools [103]. 
In the last section, we compute, as an application, the one-loop QCD corrections to the process 
pp — ► ZZZ at the LHC, also showing distributions for physically interesting quantities. 

7.2 The OPP method 

The starting point of the OPP reduction method is the general expression for the integrand of a generic 
m-point one-loop (sub-)amplitude 

= n rf {q \ ' A = (?+Pi) 2 -™i, p Q + . (56) 

In the previous equation, we use a bar to denote objects living in n = 4+e dimensions, and q 2 = q 2 +q 2 , 
where q 2 is e-dimensional and (q ■ q) = 0. N(q) is the 4-dimensional part of the numerator function 
of the amplitude. If needed, the e-dimensional part of the numerator should be treated separately, as 
explained later. N(q) depends on the 4-dimensional denominators A = (q + pi) 2 — mf as follows 

m—1 m—1 

N(q) = ^2 d(i iii 2 k) + d(q; i ii«2*3) A 
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%0<i\ i^io,ii 
m—1 m—1 

+ [a(*o) + 5(g;«o)] JJ A 

to i^io 
m—1 

+ P(q) J] A • (57) 

i 

Inserted back in Eq. (56), this expression simply states the multi-pole nature of any m-point one-loop 
amplitude, that, clearly, contains a pole for any propagator in the loop, thus one has terms ranging from 
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1 to m poles. Notice that the term with no poles, namely that one proportional to P(q) is polynomial 
and vanishes upon integration in dimensional regularization; therefore does not contribute to the ampli- 
tude, as it should be. The coefficients of the poles can be further split in two pieces. A piece that still 
depend on q (the terms d, c, b, a), that vanishes upon integration, and a piece that do not depend on q 
(the terms d, c, b, a). Such a separation is always possible and the latter set of coefficients is immedi- 
ately interpretable as the ensemble of the coefficients of all possible 4, 3, 2, 1 -point one-loop functions 
contributing to the amplitude. 

Once Eq. (57) is established, the task of computing the one-loop amplitude is then reduced to 
the algebraical problem of fitting the coefficients d, c, b, a by evaluating the function N(q) a sufficient 
number of times, at different values of q, and then inverting the system. That can be achieved quite 
efficiently by singling out particular choices of q such that, systematically, 4, 3, 2 or 1 among all possible 
denominators A vanishes. Then the system of equations is solved iteratively. First one determines all 
possible 4-point functions, then the 3-point functions and so on. For example, calling q^ the 2 (in general 
complex) solutions for which 

D = D x = A = D 3 = , (58) 

(there are 2 solutions because of the quadratic nature of the propagators) and since the functional form 
of d(q; 0123) is known, one directly finds the coefficient of the box diagram containing the above 4 
denominators through the two simple equations 

N(q±) = [d(0123) + d(^;0123)] [J A(^) • (59) 

#0,1,2,3 

This algorithm also works in the case of complex denominators, namely with complex masses. Notice 
that the described procedure can be performed at the amplitude level. One does not need to repeat the 
work for all Feynman diagrams, provided their sum is known: we just suppose to be able to compute 
N{q) numerically. 

The described procedure works in 4 dimensions. However, even when starting from a perfectly 
finite tensor integral, the tensor reduction may eventually lead to integrals that need to be regularized 
(we use dimensional regularization). Such tensors are finite, but tensor reduction iteratively leads to rank 
m m-point tensors with 1 < m < 5, that are ultraviolet divergent when m < 4. For this reason, we 
introduced, in Eq. (56), the d-dimensional denominators A, that differs by an amount q 2 from their 
4-dimensional counterparts 

Di = Di + q 2 . (60) 

The result of this is a mismatch in the cancellation of the d-dimensional denominators of Eq. (56) with 
the 4-dimensional ones of Eq. (57). The rational part of the amplitude, called R\ [170], comes from 
such a lack of cancellation. A different source of Rational Terms, called R 2 , can also be generated from 
the e-dimensional part of N(q) (that is missing in Eq. (56)). For the time being, it should be added by 
hand by looking at the analytical structure of the Feynman Diagrams of via a dedicated set of Feynman 
Rules. Examples on how to compute R 2 are reported in [170] and [171, 172]. The Rational Terms R\ 
are generated by the following extra integrals, introduced in [95,96] 

[ <^4ft 

J DiDj 

f d n~ q f 

J DiDjD k 

The coefficients of the above integrals can be computed by looking at the implicit mass dependence 
(namely reconstructing the q 2 dependence) in the coefficients d, c, b of the one-loop functions, once q 2 is 
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reintroduced through the mass shift m 2 — > m| — q 2 . One gets 

b{ij;q 2 ) = b{ij)+q 2 b {2) (ij), c{ijk; q 2 ) = c{ijk) + <fc (2) (ijk) 



(62) 



Furthermore, by defining 

m— 1 m—1 

V (m \q,q 2 )= [d(iohi2i3-,f) + d(q-,iohi 2 i3-,q 2 )] ]J A, (63) 

io<u<«2<«3 i^ioM&te 

the following expansion holds 

m 

vW(q,cf) = Y,Q {2j ~ 4) d {2j - 4 Hq), (64) 
where the last coefficient is independent on g 

d(2m-4) (?) =d (2m-4)_ (65) 

In practice, once the 4-dimensional coefficients have been determined, one can redo the fits for different 
values of q 2 , in order to determine b^ 2 \ij), c^ 2 \ijk) and S 2m ~ 4 \ Such three quantities are the coef- 
ficients of the three extra scalar integrals listed in Eq. (61), respectively. Therefore, the OPP method 
allows an easy and purely numerical computation of the Rational Terms of type Ri. 



7.3 CutTools and the problem of the Numerical Inaccuracies 

A FORTRAN 90 program (CutTools) implementing the OPP method can be found in [103], to which 
we refer for more details. We just mention that the only information needed by the code is the number 
and type of contributing propagators and the numerator function N(q) (and its maximum rank). A 
particularly interesting feature of the OPP technique, also implemented in CutTools, is that it allows 
a natural numerical check of the accuracy of the whole procedure. Given the paramount importance of 
this issue in practical calculations, we describe it here in some detail. 

During the fitting procedure to determine the coefficients, numerical inaccuracies may occur due 

to 

1) appearance of Gram determinants in the solutions for which 4, 3, 2 or 1 denominators vanish; 

2) vanishing of some of the remaining denominators, when computed at a given solution; 

3) instabilities occurring when solving systems of linear equations; 

In principle, each of these three sources of instabilities can be cured by performing a proper expansion 
around the problematic (i.e. exceptional) Phase-Space point. However, this often results in a huge 
amount of work that, in addition, spoils the generality of the algorithm. Furthermore, one is anyway 
left with the problem of choosing a separation criterion to identify the region where applying the proper 
expansion rather than the general algorithm. 

The solution implemented in CutTools is, instead, of a purely numerical nature and relies on a 
unique feature of the OPP method: the fact that the reduction is performed at the integral level. In detail, 
the OPP reduction is obtained when, as in Eq. (57), the numerator function N(q) is rewritten in terms of 
denominators. Therefore N(q) computed for some arbitrary value of q by using the 1. h. s. of Eq. (57) 
should always be numerically equal to the result obtained by using the expansion in the r. h. s. This is 
a very stringent test that is applied in CutTools for any Phase-Space point. When, in an exceptional 
Phase-Space point, these two numbers differ more than a user defined quantity, the coefficients of the 
loop functions for that particular point are recomputed by using multi-precision routines (with up to 
2000 digits) contained in CutTools [173, 174]. The only price to be payed by the user is writing, 
beside the normal ones (namely written in double-precision), a multi-precision version of the routines 
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computing N(q). The described procedure ensures that the coefficients of the scalar loop functions are 
computed with a precision defined by the user. Finally, one should mention that, usually, only very few 
points are potentially dangerous, namely exceptional, so that a limited fraction of additional CPU time is 
used to cure the numerical instabilities, therefore compensating the fact that the multi-precision routines 
are by far much slower than the normal ones. This procedure has been shown to work rather well in 
practice, as we shall see in the next section. 

7.4 pp — >• ZZZ at one-loop 

The calculation is composed of two parts: the evaluation of virtual corrections, namely one-loop con- 
tributions obtained by adding a virtual particle to the tree-order diagrams, and corrections from the real 
emission of one additional massless particle from initial and final states, which is necessary in order 
to control and cancel infrared singularities. The virtual corrections are computed using the OPP reduc- 
tion [95,96]. In particular, we make use of CutTool s [103]. Concerning the contributions coming from 
real emission we used the dipole subtraction method [175] to isolate the soft and collinear divergences 
and checked the results using the phase space slicing method [176] with soft and collinear cutoffs, as 
outlined in [177]. 

These results have also been recently presented, following a very different approach, by Lazopou- 
los et al in Ref. [8]. A more complete study, that will also include the case of W + W~ Z, W^ZZ, and 
VF + W~VF ± production, will be presented in a forthcoming publication [178]. 

Let us begin with the evaluation of the virtual QCD collections to the process qq — ► ZZZ. We 
consider the process 

?(pi) + 9(P2) — Z( P3 ) + Z(p 4 ) + Z( P5 ) (66) 
All momenta are chosen to be incoming, such that ^ pi = 0. 

Pi — - — r P$ 



P2 



Pi 



P3 



Fig. 10: Tree-level structure contributing to qq — > ZZZ. 

At the tree-level, there are six contributions to this process, obtained by the diagram illustrated in 
Fig. 10 by permuting the final legs in all possible ways. One-loop corrections are obtained by adding 
a virtual gluon to the tree-level structures, as depicted in Fig. 11. Each of the eight diagram of Fig. 11 
should be evaluated for six permutations of the final particles: overall this calculation involves the reduc- 
tion of 48 diagrams. 

We perform a reduction to scalar integrals using the OPP reduction method [95,96]. As described 
in Section 7.2, we need to provide the numerical value of the numerator of the integrand in the loop 
integrals. The numerator function N(q) can be expressed in terms of 4-dimensional denominators Di 
according to the decomposition of Eq. (57). For the particular case of five denominators, that is the 
relevant case for the process studied in this paper, we have m = 5 and the indices range from to 4. 
Next, simply by evaluating the numerator function N(q) for a given set of values of q, we can extract all 
the coefficients in Eq. (57). 

The coefficients determined in this manner should be multiplied by the corresponding scalar in- 
tegrals. Since, in the process that we are studying, no q-dependent massive propagator appears, we 
will only need massless scalar integrals. They are computed using the package OneLOop written by 
A. van Hameren [122]. 
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Fig. 11: Diagrams contributing to virtual QCD corrections to qq — > ZZZ 

As an example, let us consider the pentagon diagram (the last diagram of Fig. 11). In our notation, 
the integrand will read 

A 5 (q\ = ^(g) (67) 

[q 2 ] [(q + pi) 2 } [{q + Pi + p 5 ) 2 ] [(<? ~P2~ P3) 2 ] [{q ~ P2) 2 ] 

with 

N 5 (g) = {u( P2 ) l a P iq - P2 ) V 3 Z P iq - P2 - P3) Vi P {q+Pl+P&} V£ P iq+Pl) 1° u(pi)} (68) 

The function P(q) is the numerator of the quark propagator 

p {q) =q'+m, 

while Vf = V z • 6i , namely the contraction between he polarization vector of the i-th Z boson t\ and 
the 7-matrix in the vertex Zqq 

V» = ^ (j9ju- + 9fU+) (69) 
where ^ 

4 = ~Qf , aj = ~ S Qf , ^± = (l±7 5 )/2. (70) 

For any fixed value qo of integration momentum, and for a given phase-space point, A^s(go) is simply the 
trace of a string of known matrices. After choosing a representation for Dirac matrices and spinors, we 
evaluate N(q) by performing a naive matrix multiplication. By providing this input to the reduction algo- 
rithm, we can compute all the coefficients of the scalar integrals (in other words, the "cut-constructible" 
part of the calculation). 

The last step is the calculation of Rational Terms. As explained in Section 7.2, part of this con- 
tribution, that we call R\, is automatically included by the to the reduction algorithm. The second term 
R2, coming from the e-dimensional part of N^(q), has been added by hand by looking at the Feynman 
Diagram and turns out to be proportional to the tree-order amplitude. 

In the same fashion, we can repeat the calculation for the other seven diagrams. However, our 
method allows for a further simplification: for each fixed permutation of the final legs, only the q- 
dependent denominators of Eq. (67) will appear in the remaining diagrams. Therefore, we can combine 
all diagrams in a single numerator function and perform the reduction directly for the sum of such dia- 
grams, allowing for a one-shot evaluation of the resulting scalar coefficients. 

We checked that our results, both for poles and finite parts, agree with the results obtained by the 
authors of Ref. [8]. 

In what concerns the real emission, we only have to deal with initial state singularities, where we 
distinguish qq and qg initial states. For the qg initial state, no soft singularity is present because the 
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corresponding tree-level contribution vanishes. We recall that the structure of the NLO partonic cross 
sections is as follows: 
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(71) 



vvv 



VVVg 



where da B , da v , da c , da R , da A are respectively the Born cross section, the virtual, virtual counterterm, 
real and real-subtraction cross sections. For the qq initial state two dipoles are needed as subtraction 
terms. If p§ is the momentum which can become soft or collinear, the dipole term for gluon emission off 
the quark is given by 

8Tra s CF {1 + x 2 ' 



P<J196,<?2 



2x pi ■ pq \ 1 — x J qq 

Pl-P2-P2-P6~Pl- P6 
Pi -P2 



\MU{P})\' 



(72) 



where the {p} are. redefined momenta, {pj} = {pie, P2-, Pz-, Pa, Ps}> which are again on-shell and go to 
{pi, . . . ,p5} in the singular limit, e.g. pie = xp\. The regularised real emission part then reads 



j R _ j A 
a(J qq u<J qq 



1 1 1 



Cf \M^({ Pj })\ 1 - V qi96 ' q2 - V« 296 > qi 



'qqWVj 



VVVg , 



where the factor 1/6 accounts for the three identical bosons in the final state. More details can be found 
in [175,178]. 

The hadronic differential cross section with hadron momenta Pi and Pi is the sum over all partonic 
initial states convoluted with the partem distribution functions 



da(P 1 ,P 2 ) = J dz 1 dz 2 fa(z 1 ,HF)Mz2,PF)da ab (z 1 P 1 ,z 2 P2) 

ab 



(73) 



where the sum runs over the partonic configurations qq, qq, gq, qg, gq, qg. 



7. 4. 1 Numerical results 

As an explicit example we present the numerical results for the case uu — > ZZZ for y/s = 14 TeV 
and using CTEQ6L1 [28]. The tree-order cross section has been evaluated using the HELAC event 
generator [179-181]. In the following table the results in fb are presented for the tree-order cross section 
<7o, the ratio of the virtual to the tree-level cross section, and the real contribution, combining 5— and 
6— point contributions, as described above, for all channels, i.e., uu,ug,gu, for different values of the 
factorization(renormalization) scale (fi = fip = hr). 



scale 


O"0 


ay/a 


&R 


&NLO 


A* 


= M Z 


1.481(5) 


0.536(1) 


0.238(2) 


2.512(2) 


A* 


= 2M Z 


1.487(5) 


0.481(1) 


0.232(2) 


2.434(2) 


A* 


= 3M Z 


1.477(5) 


0.452(1) 


0.232(2) 


2.376(2) 


A* 


= 4M Z 


1.479(5) 


0.436(1) 


0.232(2) 


2.355(2) 


A» 


= 5M Z 


1.479(5) 


0.424(1) 


0.237(2) 


2.343(2) 



As it is evident from these results, the K— factor is quite sizable (1.58 — 1.69), whereas the 
dependence on the scale \i is for both cases quite weak, due mainly to the electroweak character of the 
process. 
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7.5 Conclusions 

We presented a new method for NLO processes (OPP), in which the reduction to known integrals is 
performed at the integrand level. The method has been successfully tested in a number of applications, 
the latest being the production of three Z bosons at the LHC. 

The efficiency of the method is quite good. It can be further improved if the numerical evaluation 
of the integrand in the one-loop amplitude, by means of recursion relations, without relying on Feynman 
diagrams, is developed [67]. 

In general, the speed, the precision and the simplicity of the OPP method, make it a very good 
candidate for the construction of a universal NLO calculator/event-generator. 

Part II 

IMPROVEMENTS ON STANDARD 
TECHNIQUES 

8. GOLEM: A SEMI-NUMERICAL APPROACH TO ONE-LOOP AMPLITUDES 14 
8.1 Introduction 

The first collision data from the Large Hadron Collider (LHC) at CERN are expected in a couple of 
months, giving us the opportunity to explore unprecedented energies and luminosities. However, in 
order that a discovery of New Physics can be claimed, it is of crucial importance to have the Standard 
Model physics under control. This includes e.g. understanding of the detectors, the underlying event, 
the luminosity determination, the jet energy scale [182]. For most of these issues, an inteiplay between 
measurements and precise theory predictions is mandatory. In a hadron collider environment, multi- 
particle/jet final states will be produced in abundance. Therefore considerable effort needs to be spent 
to make predictions for multi-particle processes beyond the leading order. While the calculation of one- 
loop five-point amplitudes can be considered as the state of the art at the moment, the first complete cross 
section for six-point processes at hadron colliders still awaits its completion. Many different approaches 
to multi-particle production have been developed in the last few years, most of them being described in 
these proceedings. For other reviews and very recent developments, see e.g. [103, 132, 183]. 

Here we will focus on a method implemented in the program GOLEM (General One-Loop Eval- 
uated of Matrix elements), which is based on a semi-numerical evaluation of building blocks stemming 
from the reduction of one-loop Feynman diagrams [184]. The main features of the formalism are the 
following: 

• It is valid for massive and massless particles 

• For N > 5 external legs, the reduction of rank R iV-point integrals is done algebraically, reducing 
the rank and the number of propagators at the same time in each reduction step. For N < 5 
we worked out form factor representations which allow to avoid inverse Gram determinants in 
exceptional kinematic regions. 

• The infrared divergences are easily extracted analytically in terms of triangles. 

• The rational parts of the amplitudes are obtained as by-products and can be projected out. 

• The program has an analytic and a numerical branch: it can perform a complete reduction to 
scalar integrals, represented in terms of analytic functions: such a complete reduction introduces 
inverse Gram determinants, but this branch can be chosen safely in phase space regions where the 
Gram determinants are sufficiently large (which is the bulk of the phase space). As the evaluation 

"Contributed by: C. Bernicot, T. Binoth, J.-Ph. Guillet, G. Heinrich, E. Pilon, T. Reiter 
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of analytic functions is fast, this speeds up the program considerably as compared to a purely 
numerical approach. Near exceptional phase space points, the program allows to stop the reduction 
before dangerous denominators are produced. The building blocks to evaluate in this case are 
finite three- and four-point functions with Feynman parameters in the numerator. As a brute- 
force numerical evaluation of the four-point functions is rather slow, we have worked out one- 
dimensional integral representations, whose numerical evaluation is extremely fast. Details will 
be given in the following section. 

We have implemented the reduction in algebraic manipulation programs and have obtained fully an- 
alytical results for several amplitudes using these methods [2,94, 185-188]. Without having efficient 
and automated simplification methods to reduce the size of big analytic expressions, the fully analytic 
approach based on form factors suffers from factorial complexity and therefore does not seem to be 
appropriate for 6-point processes. The semi-numerical reduction is preferable in this case. For the cal- 
culation of the rational terms alone the situation is different, as the form factor representations simplify 
considerably when restricted to terms which can generate rational parts [93]. 

8.2 Results 

Below we will describe applications of our method to one-loop six-point amplitudes and explain in detail 
certain features which guarantee a fast and numerically robust evaluation in all phase space regions. 

8.2.1 The GOLEM numerical library 

In the GOLEM library, the strategy is to evaluate numerically higher dimensional three- and four-point 
functions in phase space regions where numerical instabilities arise due to spurious singularities. To be 
specific, these integrals are six- and eight-dimensional four-point functions I® +2 , and four- and 

six-dimensional three-point functions 1$ , I% +2 , with or without Feynman parameters in the numerator. 
While the triangles are two-dimensional integrals in Feynman parameter space, the boxes a priori involve 
integration over three Feynman parameters. As numerical integrations in multi-dimensional parameter 
space are rather slow, we worked out one-dimensional integral representations for these integrals, whose 
evaluation is both fast and precise. In [184, 189] we have already presented other methods for the numer- 
ical evaluation of Feynman parameter integrals, but the one-dimensional representations discussed here 
are preferable, as they are much faster. 

As an example, let us consider the case where two massive particles scatter into two light particles 
via a fermion loop. The two ingoing particles have a small velocity. In this kinematic region, the Gram 
determinant is small. In this case, we have to evaluate four-point functions with two adjacent massive 
legs, and with Feynman parameters in the numerator. In fig. 13 we plot the six-dimensional four-point 
function with two adjacent massive legs, /| ad j, against the absolute value of the coefficient B which is 
proportional to the ratio det(G)/det(S), for a trajectory of points with 10~ 15 < \B\ < 1CT 3 . 

In the GOLEM library, there is a cut c which allows to split the phase space regions where the 
four-point function is evaluated analytically from those where it is evaluated numerically. The larger the 
cut, the longer the evaluation takes, as more calls of the numerical integration routine are made. On the 
other hand, if the cut is too small, the analytical evaluation causes a loss of precision of several digits. 

As an illustration, we compute T\{z\z^) and we plot the real and imaginary parts for different 
values of the cut c: c = 1CT 1 (Fig. 13), c = 1CT 3 (Fig. 14) and c = 1(T 5 (Fig. 15). 

In the case at hand, the CPU time does not vary very much with the cut, the evaluation time ranges 
from 0.14 s (on an Intel Pentium M 1.3 GHz) for c = 1CT 1 to 0.10 s for c = 10~ 5 . However, this 
statement is hard to generalise to all possible situations occurring in a calculation of a complex multi-leg 
amplitude. In any case, the cut c allows to adjust the trade-off between speed and precision. 
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8.2.2 The uu — ► ddss amplitude 

With our method we calculated the one-loop six-quark amplitude 

A(u(pi, \i),u(p 2 , A 2 ) -> d(p3, A 3 ),(J(p4, A 4 ),s(p5, A 5 ),s(p 6 , A 6 )) (74) 

in massless QCD. The calculation has been carried out using spinor helicity amplitudes in the 't Hooft- 
Veltman scheme. We have chosen a convenient colour basis, which allows to split the amplitude as 
follows 

6 

^^C^( m ,...,p 6 ), (75) 

A i=l 

where A c are the helicity and colour subamplitudes. In particular we chose the colour structures 

(C\ C 2 , C 3 , C 4 , C 5 , C 6 ) = (S%6*S%, 8%8«S%, 8*8*6%, 5%6%5%, ( 76 ) 

In our notation A is the vector (Ai, . . . , and Aj = ±1 is the helicity of the particle with momentum 
Pj of which the colour index is Cj. In the six-quark amplitude one can identify two independent helicites 
A a = (+, +, +, +, +, +) and X b = (+, +, +, +, — , — ); all other helicities are either identically zero or 
related to A a or X b by parity invariance, which is exploited in our calculation. 

We generated the Feynman diagrams for this process with QGraf [165] and reduced the tensor 
integrals using FORM [166, 190] to form factors as defined in [184]. We deal with the spinor algebra by 
completing spinor lines to traces, e.g. for an arbitrary product F of Dirac matrices we use 

<p ' +|r i^) = %4s> ,r{(i+75)i<j ^' r} ' <77) 

With the help of FORM and Java code the expressions for the diagrams are transformed into a Fortran90 
program. The Golem9 library is used for the numerical evaluation of the form factors. In this approach 
we found it advantageous to treat the spinor traces numerically as well, in order to keep the expressions 
more compact. 

The code returns the subamplitudes in the form 

AKpi, ■ ■ • ,P6) = ^2~ s (A + f + <?+ 0(e)) (78) 

for each of the six colour structures and for all non-zero helicities, where A, B and C are complex 
coefficients. As an example we plot in Figure 16 the quantity s\ A^\a~^ for one colour structure C 1 and 
the two helicity configurations A a and X b . The initial state momenta are chosen to be along the z-axis 
while the final state momenta have been rotated about the y-axis by an angle 9. For 9 = the momenta 
are chosen as in Ref. [191]: 

p 3 = (33.5,15.9,25.0) 
p 4 = (-12.5,15.3,0.3) 
p 5 = (-10.0,-18.0,-3.3) 

p 6 = (-11.0,-13.2,-22.0) (79) 

In the chosen units the renormalisation scale is /x = 1. The amplitude has been evaluated at 50 successive 
points between 9 = and 9 = 2ir (9 = 0, 0.126, 0.252, . . .), which took 2.4 seconds per point and 
helicity on an Intel Pentium 4 CPU (3.2 GHz). Table 3 shows the numerical values of all coefficients for 
the point 9 = 0. 
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_c A 

1 -0.0029670 - 0.0036065i 

2 0.0042784 + 0.0049474z 

3 -0.0123663 - 0.0186981i 

4 0.0051836 + 0.0066459i 

5 -0.0143367 - 0.0137603i 

6 0.0083400 + 0.0100456i 



B 

0.0203701 + 0.0281510i 

-0.0191448 - 0.0420120i 
0.1171088 + 0.1401148i 

-0.0462621 - 0.0477458i 
0.1282264 + 0.1049820i 

-0.0745825 - 0.0730179i 



C 

-0.0659100 - 0.1057940i 

0.0338141 + 0.1820798i 
-0.4902357 - 0.4754639i 

0.1803702 + 0.1706208i 
-0.5199953 - 0.3972433i 

0.2929410 + 0.2459317f 



Table 3: Six-quark amplitude. Numerical values of the virtual part A* {{pj}j=i...6)a s 3 for the kinematics given in the text 
and 6 = 0. 




Fig. 12: Scale-dependence of the cross section onlo (solid line) compared with the tree-level cross section cto (dashed line). 
The scale is reported in the plot in units of Mz, from /j, = Mz to /j, — 5Mz. 
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Fig. 13: The six-dimensional four-point function with three Feynman parameters in the numerator, if (212I), with two adjacent 
massive legs and the cut c — 10 _ . 
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Fig. 14: The six-dimensional four-point function with three Feynman parameters in the numerator, I^ziz^), with two adjacent 
massive legs and the cut c = 10~ 3 . 
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Fig. 1 5 : The six-dimensional four-point function with three Feynman parameters in the numerator, I±(z\Z2), with two adjacent 
massive legs and the cut c = 1CP J . 
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e 

Fig. 16: The six-quark amplitude. The finite parts of the Laurent expansion in e of s\Ai \aj 3 (solid) and s\Ai \aj 3 (dashed) 
are plotted for a kinematic point defined in the text, where the final state momenta have been rotated about the y-axis by an 
angle 6. 
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9. ISSUES WITH THE LANDAU SINGULARITIES 15 
9.1 Introduction 

Cross sections involving a large number of external particles can contain numerical instabilities which 
must be carefully located and controlled. At tree-level one can mention integration over a i-channel 
pole if the integration variables are not properly chosen. The crossing of a resonance might also be 
problematic. Beside these physical situations there might be fake singularities specific to the way one 
has set up the amplitude; one example is the singularity brought about by an unlucky choice of a reference 
vector at the helicity amplitude level. These problems are exacerbated at the loop level since the loop 
integrals can also develop singularities. A prominent example is the occurrence of vanishing inverse 
Gram determinants: see for example the contribution of Denner and Dittmaier. The latter is a fake 
singularity that can be met for some special, and simple, kinematical conditions on the phase space 
of the external particles having to do with how one has chosen one's (independent) basis for the loop 
integrals and how one has subsequently expressed the other loop integrals in this basis. Loop integrals 
can also have true singularities that have an underlying physical origin. They depend on the dynamics 
of the problem. Thresholds are one example, though harmless and trivial to locate. These types of 
singularities belong to the general class of Landau singularities. The physical singularity can be revealed 
by studying the analytic properties of the scalar integral. Here we study the case of one-loop integrals. 
In particular we will review how the conditions for having such singularities can be derived, especially 
in a format that is conducive to an easy implementation in a computer code. When such a singularity 
is present it is important to inquire whether this singularity is integrable or not. We rederive here the 
singular part. We then consider two specific complementary examples taken from the recent literature. 
The first one, the electroweak corrections to pp — > bbH, reveals a Landau singularity having to do with 
massive, indeed unstable, particles in the loop. In this case the singularity is smoothed out by the width 
of the unstable particles. The second is the 6-photon amplitude which involves massless states, both 
internally and externally. In this case the Landau determinant is a quadratic function whose square root 
is proportional to the Gram determinant. 



9.2 Conditions for a Landau singularity and the nature of the singularity 

Consider the one-loop process Fi(pi) + i 7 ^^) + 
. . . + Fn(pn) — ► 0, where F{ stands for either a 
scalar, fermion or vector field with momentum pi 
as in the figure opposite. The internal momentum 
for each propagator is q\ with i = 1, . . . N. Each 
momentum qi is associated with one Feynman pa- 
rameter Xi respectively. The scalar loop integral 
reads 



rpN _ 

I o = 



Di 



d D c 



1 



qi-m^ + ie, qi = q + n, 



i=i 



1,...,N, 



(80) 
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The Feynman parameter representation reads 



poo N p jD i 

T N = T(N) / dx 1 ---dx N 8{Y J x i -l) / * - - — ? . (81) 

JO ~{ J (2tt) u i{xiDi+X 2 D2^ X N D N )" 

Because of the Dirac delta function, the integration boundary in the Feynman parameter space are X{ = 0, 
i = 1, . . . , N. Thus the only important condition on xi is that they are not negative. The singularities are 
given by the Landau conditions [123, 192] 

Vi Xi(qf - ml) = 0, 

Li=l x iQi = °- 

If eq. (82) has a solution X{ > for every i 6 {1, . . . , N}, i.e. all particles in the loop are simultaneously 
on-shell, then the integral has a leading Landau singularity (LLS). If a solution exists but with some 
Xi = while the other Xj's are positive, the Landau condition corresponds to a lower-order Landau 
singularity (LOLS). 

By introducing the matrix Q, under the condition qf = mf, 



it %'■ 

Qij = 2 qi . qj = m 2 + m) - ( ft - qj f = mf + m) - (n - r^) 2 ; i,j G {1, 2, . . . , M}, (83) 
the conditions to have a Landau singularity in the physical region are 



det(Q) = 0, 

Xi > 0, i = 1, ... ,M. 



(84) 



For M = N one has a leading singularity, otherwise if M < N this is a subleading singularity. If some 
internal (external) particles are massless, as in the case of six-photon scattering, then some Qij are zero, 
and the above conditions can be easily checked. However, if the internal particles are massive then it is 
difficult to check these conditions explicitly, especially if M is large. In this case, we can rewrite the 
above conditions as follows 

det(Q) = 0, 

Xj = det(Q jM )/ det(Q M u) > 0, j = 1, . . . , M - 1, 

where Qij is obtained from Q by discarding row i and column j from Q. Note that det(QjUAf) = 
d[det(Q)]/dQMM- If det(QMAf) = then the second condition in (85) becomes det(Q_jjw) — 
with j = 1, . . . , M — 1. There may be cases, as we will encounter in section 9.4, where the Landau 
determinant det(Q) has a quadratic form. These special situations have to be handled with care. 

The existence of a Landau singularity corresponds to an eigenvector of Q with zero eigenvalue. In 
general, Q has N real eigenvalues Ai, . . . , Ajv- Consider the case where Q has only one (non-degenerate) 
very small eigenvalue Ajv -C 1. To leading order 

Xn = — , a\ = A1A2 • • • Ajv-i ^ 0, a = det(Q). (86) 
ai 

With V = {x®,X2, ■ ■ ■ ,x° N } the eigenvector corresponding to Ajy, we define v 2 = V ■ V. We will 
assume that Ai > for % = 1, . . . , K and Xj < for j = K + 1, . . . , JV - 1 with < K < N - 1. It 
can then be shown that in D dimensions, 

(.^^(A^-p^ ^- -i)/2r ((jV _ g + 1 )/2) 

2 3D / 2 - J V(-l) Ar-i ^ 1 ai (A^ 2 -ie)^ D+1 )/ 2 ' 
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This result holds provided a\ / and N - D + 1 > 0. For the box, N = 4, D = 4, o 
we get 

e i7r(3-/0/2 



(r 4 : 



4V(-l) 3 -^det(Q 4 ) -ie' 
This shows that (Tq)cHv is integrable but its square is not. 
In the case N = 3 (the triangle), D = 4, it is possible to derive 

_ e i7r(2-K)/2 u 

ln(A3i; — ie) 



and a± ^ 
(88) 



T 3 



8V(-l) 2 -^AiA 2 



(89) 



Tq and its square are therefore integrable. 



9.3 



99 



bbH 



The first example we study is the electroweak corrections to pp — ► [193] where the one-loop 
amplitude squared, which is all that remains in the limit of vanishing bottom Yukawa coupling, develops 
a Landau singularity which represents the rescattering of the top pair and their decay into a W pair that 
produces the Higgs through WW fusion. As we will see, in this example, introducing the width of the 
internal top and W particles smoothes the singularity. There is a leading Landau singularity present in 
the box diagram shown in Fig. 17 that occurs for some specific values of the kinematic variables. 




00000"0"0"0TJ 




Fig. 17: A box diagram contributing to gg — > bbH that can develop a Landau singularity for AIh > 2Mw and yfs > 2mt 
i.e. all the four particles in the loop can be simultaneously on-shell. 



With g{pi)+g(p 2 ) b( P3 )+b(p 4 )+H(p 5 
and the on-shell conditions p\ = p\ = 0, p\ 
box integral is a function of two variables s\p, 

T 4 (si, s 2 ) = A)(Mjr, 0, s, 0, si, a 2 , m 2 



s = (pi+p 2 ) 2 , si = (P3+P5) 2 , s 2 = (P4+P5) 2 , 



0, = Mfj, fixing s and Mh, the scalar 



(90) 



The kinematically allowed region is 



Mjj < si < s, Mfj— <s 2 < Mir + s- si 
si 



(91) 



The reduced matrix, which is equivalent in this case to the Q matrix for studying the Landau 

singularity, is given by 



ij 



2M^- 



M% 



si 



2M w m t 
\ 2 M w m t 



2 M w m t 
2M w m t 



2 M w m t 

2M w m t 
1 

2m, 2 -s 
2 mi 



2M w m t 

2M w m t 
2fiij-s 
~2^f 

1 



(1) 



Q 



2mjrrii 



(92) 
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The singularity corresponds to det(5^ 4 ^) = 0. The determinant is a quadratic function of si,s% when 
s and all internal masses are fixed. The Landau determinant, the real and imaginary parts of Tq are 
displayed in Fig. 18 for yfi = 353 GeV, M H = 165 GeV, m t = 174 GeV, M w = 80.3766 GeV. 
We clearly see that the Landau determinant vanishes inside the phase space and leads to regions of 
instability exhibiting leading and lower-order Landau singularities in the real and imaginary parts of 
the scalar integral. To investigate the structure of the singularities in more detail let us fix ^fs\ = 





Fig. 18: The Landau determinant as a function of Si and S2 (upper figure). The real and imaginary parts of Dq as a function 
of si and S2. 



^2{m% + Myy) ~ 271.06 GeV, so that the properties are studied for the single variable «2- The results 
are shown in Fig. 19. 

From Fig. 19 we see that there are four discontinuities in the function representing the real part 
of the scalar integral in the variable y/~S2. As S2 increases we first encounter a discontinuity at the 
normal threshold J~S2 = m t + My/ = 254.38 GeV. This corresponds to the solution (for the Feynman 
parameters) = and X2,4 > of the Landau equations. The second discontinuity occurs at the 
anomalous threshold ^/~S2 = 257.09 GeV of a reduced triangle diagram. This corresponds to the solution 
X3 = and > of the Landau equations. The condition of vanishing determinant det(5 < 3) = for 
this triangle has two solutions 

S2 = 2^2 - (m£K + M&) T M Hy /M>-4Mfa(m? t - M^)j (93) 
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Fig. 19: The imaginary, real parts of Do and the Landau determinant as functions of S2- 



which gives = 257.09 GeV(inside of phase space) and 297.86 GeV (outside of phase space). We 
can also check that the former value satisfies the sign condition in (85) while the latter does not. Note 
that one of the conditions for this anomalous threshold to occur in the physical region is Mh > 2M\y, 
see Eq. (93). The same phenomenon happens for the third discontinuity at ^fs2 = 259.58 GeV 
which corresponds to the anomalous threshold of the reduced three point function obtained from the 
box diagram by contracting to a point the x\ line. The last singular discontinuity is the leading Lan- 
dau singularity. The condition det(S , 4) = for the box has two solutions which numerically corre- 
spond to yf~S2 = 263.88 GeV or ^/s2 = 279.18 GeV. Both values are inside the phase space, see 
Fig. 19. However after inspection of the corresponding sign condition only = 263.88 GeV (with 
xi « 0.533186, x 2 « 0.748618, x 3 « 0.774941) qualifies as a Landau singularity. ^ = 279.18 GeV 
has x\ « —0.742921, x 2 ~ —0.748618, x$ « 1.06537. The nature of the leading Landau singularity in 
Fig. 19 can be extracted by using the general formula (88). With the input parameters given above, the 
Landau matrix has only one positive eigenvalue at the leading singular point, i.e. K = 1. The leading 
singularity behaves as 



16M^mfVdet(£ 4 ) -is 

When approaching the singularity from the left, det(54) > 0, the real part turns singular. When we 
cross the leading singularity from the right, det(54) < 0, the imaginary part of the singularity switches 
on, while the real part vanishes. In this example, both the real and imaginary parts are singular because 
det(54) changes sign when the leading singular point is crossed. 

The instabilities of the integral and the singularities are due to the unstable internal particles. The 
problem can be remedied by introducing the finite width of the W and top. As seen from Fig. 20, 
introducing the finite width effect in the scalar box gives a smooth behaviour. 

9.4 The six photon amplitude 

The second example concerns a case with massless internal particles involving massless external parti- 
cles: the 6-photon amplitude [194], see also [94,96, 191]. Although the scalar integrals for the 6-photon 
amplitude have a potential Landau singularity that leads to some characteristic patterns of the amplitude, 
direct calculations of the helicity amplitudes show that the singularity is tamed by the dynamics of the 
gauge interaction in a somehow unexpected way. This is welcome since we would not able, in this case, 
to revert to the trick of introducing a width for the particles. This said, introducing non-zero (internal) 
masses, as would be fit for the couplings of the massless photons, would regulate a vanishing Landau 
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Fig. 20: Effect of the combined width of the W, Yw and the top, F — t, to the real and imaginary part of the scalar function. 



determinant, but would of course still pose a considerable numerical problem if the singularity from the 
vanishing Landau determinant is not counterbalanced by the spin and gauge algebra. 

To be able to see the cancellation at the level of the amplitude is only possible if one has very 
compact analytical expressions for these amplitudes. In our investigation the expressions for the ampli- 
tudes [194] are based on the unitarity-cut methods and are made particularly simple thanks to the fact that 
the six-photon amplitude has no IR/UV divergences and no rational terms. The six-photon amplitude was 
calculated in three models: i) scalar QED, A§ calar , ii) spinor QED: j[f ermion an( j jjj) supersymmetric 
QED M = 1: A# =1 . The three amplitudes Af alar , A{ ermion and A{f =1 are in fact related through: 

A f 6 ermion = -2A s 6 calar + A% =1 (95) 



Full compact expressions for the amplitudes can be found in [194]. The potential Landau singular- 
ity in the 6-photon amplitude reveals itself in the so-called double parton scattering configuration [191], 
see Fig. 21. 




Fig. 21: Double parton scattering configuration: pi,P4 are incoming photons with pi + pi = 0, each splits into a fermion pair which 
rescatters to give photon pairs (p2 , P6 ) ; (P3 , P5 ) at very small, vanishing, transverse momentum. 



The Landau conditions read 

det(Q) = (S135S435 - S35S26) 2 -> , S35,S 2 6>0 , Si 35 ,S435 <0 (96) 

where Sijk = {pi + Pj + Pk) 2 , all the pj's are taken as incoming. Note the specific nature of det(Q) 
which has a quadratic form. This will lead to a double root (eigenvalue) at the singularity, or in other 
words the derivative of det(Q) at the singularity is also vanishing. In fact det(Q) is proportional to the 
square of the Gram determinant, det(G). To wit 

det(G) = -2si 4 (si35S 4 35 - S35S26) oc \J det(Q) (97) 
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This property is due to the presence of many zeros, both from the kinematics of the external photons and 
the masslessness of the internal lines. 

How does the singularity of the scalar integral transpire at the level of the amplitude? Let us turn to 

the NMHV ( six-photon helicity amplitude and specialise to the kinematics 16 of the Nagy 

and Soper configuration [191]. We start from a fixed point in phase space in the centre of mass frame 
Pi + Pi = with P4 along the z-axis: 

f p£ = (-33.5, -15.9, -25.0) pt = (11.0, 13.2, 22.0) 

\ ps = (12.5,-15.3,-0.3) p| = (10.0, 18.0, 3.3) 1 ; 

One can generate new configurations by rotating the final state about the y-axis by an arbitrary angle 9. 
We can then study the behaviour of the amplitude in this parameter. It is illuminating to rewrite det(Q) 
in terms of this parameter for this particular configuration: 

det(Q) = (sukf) 2 with k 2 = p\^ y + (p 35x cos 9 + p 35z sinf?) 2 (99) 

where p 35 » = p 3i + p 5i , i = x, y, z. The minimum value of k t is given by k 2 min = p\ hy . 

The behaviour of the amplitude as a function of 9 for this particular configuration is shown in 
Fig. 22. The important conclusion to draw from Fig. 22 is that the structure of the amplitude, in particular 




e e 



Fig. 22: The NMHV amplitude as a function of 9 in the Nagy-Soper configuration in the case ofQED (left) as well as the scalar and J\f = 1 
SUSY(right). In the first panel we also show the dependence ofk^ which is a good measure of dct(Q). 

the peculiar dips, is well tracked by det(Q). Indeed the dips that show in the amplitude occur exactly at 
the points where det(Q) is smallest. The dips occur at 9 ~ 2.32 and 9 ~ 2.32 + ir ~ 5.46. These values 
can be derived from Eq. (99) where k t = kt m i n - 

One can ask what would happen in a configuration where k tm i n and consequently det(Q) — ► 0? One 
can arrive at this det(Q) — ► configuration by perturbing the original kinematics in Eq. 98 

f P2 - 2 = ("33.5, -15.9 - A y , -25.0) ^ - g = (11.0, 13.2 + Aj,, 22.0) 
1 P5 ^ P5 = ( 12 -5, -15.3 + A y , -0.3) ft 2 = (10.0, 18.0 - A y , 3.3) 

The 9 modulation is unchanged, such that the dips occur at the same location in 9. However now A y can 
be chosen such that ktmin = 0. This occurs for A y = 1.05. 

Figs. 23 show how the pattern of the amplitude, as far as the dip around the singularity at 9 = 5.46 
is concerned, evolves as A y is varied from zero to 1.05 where det(Q) and k tm i n vanish. It can be seen 

16 The correspondance between the kinematical conventions of Nagy and Soper and the one used here are the following: 
Nagy and Soper consider the reaction 7 + (p4) + 7~(pi) — » 7~(— Pe) + 7 + (~ P2) + 7 + (~ P3) + 7~(— Ps) i- e - their ki's and 
ourpj are such that: ki = p4, k2 = pi, ks = — P2, &4 = — Ps, fcg = — pe and ke = — P3 so that ki + &2 = ^3 + ^4 + ^5 + fee- 
See [194] for more details. 
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that as det(Q) — > with increasing A y , the width of the dip decreases more and more so as to behave 
as a sudden jump, with the oscillation pattern disappearing completely for A y = 1.05. The numerators 
of the six-photon amplitudes, reflecting the dynamics of the gauge interaction, vanish fast enough as the 
Landau singularity is approached. Therefore the singularity seems to be dynamically regulated for the 
three cases of the scalar, the fermion and the SUSY-amplitude. 

It is also revealing to investigate how the apparent Landau singularity is approached from different 
directions by considering a two-dimensional parameterisation of det(Q) and the kinematics. 

We therefore modify the original Nagy-Soper parameterisation such as to generate a Landau sin- 
gularity and add a k t variable both along the x and y direction to follow the approach to the singularity: 

^ = (-33.5 - k tx , -15.9 - k ty , -25.0) ^ = (-12.5 + k tx , 15.3 + k ty , 22.0) 
^ = (12.5 + kt Xi -15.3 + k ty , -0.3) p£ = (33.5 - k tx , 15.9 - k ty , 3.3) 

Figs. 24 show the six-photon amplitudes ^ ermton l- N '- 1 as f unc tions of the two variables k tx and kt y . 
Up to an overall rotation, the analytic structure of these amplitudes near the Landau singularity at kt x = 
k ty = can be modelled as 

A, ~ tJ^t = ^in(2o) (102) 

K tx~*~ K ty z 

where k tx = k t cos a, kt y = h sin a, k t = (kf x + kl y ) 1 / 2 . 

The amplitudes exhibit a valley and a ridge along mutually perpendicular axes crossing each other 
at kt x = kt y = 0. The various profiles shown in Fig. 23 are nothing but cross sections at fixed k t y of 
Fig 24. In particular, the profiles for A y = 1.05 correspond to kty = 0. More generally, when both 
ktx and kty approach simultaneously, Aq remains finite: the Landau singularity of the double parton 
scattering type does not lead to a divergence ~ l/fc| as would have been naively expected from a general 
power counting argument [192, 195]. Yet the limiting value of A§ depends on the direction a along 
which the origin k t x = kt y = is approached. 

9.5 Conclusions 

We foresee that in the calculations of multi-leg one-loop processes the study of the Landau conditions will 
bring very useful, if not crucial, information. More investigations of the properties of these singularities 
need to be performed. 
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Fig. 24: The six-photon amplitude in spinor QED (top) and in J\f = 1 QED (bottom) around the Landau singularity 

10. TENSOR ONE-LOOP INTEGRALS IN EXCEPTIONAL PHASE-SPACE REGIONS 17 
10.1 Introduction 

At the LHC and ILC, many interesting processes involve more than four external particles. A thorough 
description of such processes requires the evaluation of strong and electroweak radiative corrections 
at least in next-to-leading order (NLO). The most complicated part in such calculations concerns the 
numerically stable evaluation of the one-loop tensor integrals of the virtual corrections. 

For processes with up to four external particles the classical Passarino-Veltman (PV) reduction 
[167], which recursively reduces tensor to scalar integrals, is sufficient in practically all cases. This 
scheme, however, involves Gram determinants in the denominator, which spoil the numerical stability 
if they become small. With up to four external particles this happens only near the edge of phase space 
(forward scattering, thresholds). With more than four external particles, Gram determinants also vanish 
within phase space, and methods are needed where Gram determinants can be small but still non-zero. In 
this context it should be noticed that the described problem of inverse Gram (and related) determinants 
occurs in all methods that reduce loop diagrams or amplitudes to the basis set of standard scalar integrals. 
This, in particular, also applies to unitarity-based or bootstrap approaches that work at the analytical (see 
e.g. Ref. [132] and references therein) or numerical [95,102,105,126,183] level. These methods certainly 
mitigate the problem of cancellations, but cannot avoid it completely. 



''Contributed by: A. Denner, S. Dittmaier 
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In this article we inspect two benchmark phase-space points that are inspired from our calculation 
of electroweak (EW) 0(a) corrections to e + e~ — > 4fermions [196, 197]. 18 One of the two points 
involves a small Gram determinant, the other involves both a small Gram and a small "modified Cayley 
determinant" at the same time. Although of course the real performance of proposed solutions can be 
only be found out in full applications, i.e. when integrating loop corrections to complicated processes 
over the whole phase space, a selection of such benchmark points is certainly a useful testground in the 
development of loop techniques. 

Several solutions to the problem of numerical instabilities due to inverse Gram determinants have 
been proposed in recent years, but not many of them have proven their performance in complicated 
applications yet. For references and descriptions of some methods alternative to ours, we refer to 
Refs. [200,201]. 



10.2 Tensor coefficients and their reduction 

We consistently follow the notations and conventions for scalar and tensor one-loop integrals introduced 
in Refs. [200,202]. Here we briefly repeat the conventions for 4-point integrals as required in the con- 
sidered examples. Tensor 4-point integrals of rank P are defined as 



(27T/i) 



A-D 



17T Z 



q 



Mi . . . a^P 



N k = {q+p k f -m| + i0, p = 0, 



(103) 



where D is the number of space-time dimensions and \i the reference scale of dimensionional regular- 
ization. The tensor integrals are decomposed into covariants as follows, 



3 



11,12=1 



=i 

£ jSWw + £or< +p"X + 9 Pfi Pi)Do 0il ; 



(104) 



11,12,13 = 1 



H=l 



and so on for higher rank. Up to rank 3, and only those are considered below, 4-point tensor integrals are 
UV finite. The kinematical arguments of the coefficients D. .., which comprise all scalar products piPj 
and internal masses m^, are written as 



D... = D,„(pI,(p2 -Pi) 2 ,(P3 -P2) 2 ,pi,p|,(P3 -Pi) 2 ,"io,m?,m|,m|). 



(105) 



Conventional PV reduction [167] expresses the rank-P 4-point coefficients in terms of lower-rank 4- and 
3-point coefficients. In each step P — > (P — 1) the inverse of the Gram matrix 



2pipi 2pip 2 2pip 3 
2p 2 pi 2p 2 p2 2p 2 p 3 
2p 3 p 1 2p 3 p 2 2p 3 p 3 



(106) 



occurs, which causes the above-mentioned numerical problems if the determinant \Z\ becomes small. 
The highest negative power of \Z\ occurs in the calculation of tensor coefficients Di 1 i 2 ___ without "0" 
indices, rendering them numerically the most delicate. In the following we also need the matrix 



/ 2ml 



X 



fi 
h 

V h 



h f'2 h \ 
J 



Z 



t 2 2,2 

fk=P k - m k + m . 



(107) 



18 Meanwhile the same methods have been successfully applied to NLO EW and QCD corrections to the Higgs decay H 
WW/ZZ — ► 4/ [198, 199] and to Higgs production via vector-boson fusion at the LHC [6,7]. 
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The vanishing of the modified Cayley determinant \X\ corresponds to necessary conditions for true 
(Landau) singularities in a Feynman diagram. The minors (i.e. determinants of submatrices where row 
i and column j are discarded) of the matrices Z and X, respectively, are called Zij and Xij in the 
following. 

10.3 The "DD" approach 

One-loop tensor integrals can be naturally grouped into three categories, which we have treated in com- 
pletely different ways: 

(i) For 1- and 2-point integrals of arbitrary tensor rank, numerically stable analytical expressions 
are presented in Ref. [200] (see also Ref. [167]). 

(ii) For 3- and 4-point tensor integrals, PV reduction [167] is applied for "regular" phase-space 
points where Gram determinants are not too small. For the remaining problematic cases special reduction 
techniques have been developed [200]. 

One of the techniques replaces the standard scalar integral by a specific tensor coefficient that can 
be safely evaluated numerically and reduces the remaining tensor coefficients as well as the standard 
scalar integral to the new basis integrals. In this scheme no dangerous inverse Gram determinants occur, 
but inverse modified Cayley determinants instead. We note that the procedure is related to the fully 
numerical method described in Ref. [203]. 

In a second class of techniques, the tensor coefficients are iteratively deduced up to terms that 
are systematically suppressed by small Gram or other kinematical determinants in specific kinematical 
configurations. The numerical accuracy can be systematically improved upon including higher tensor 
ranks. In our previous applications the highest relevant tensor rank was improved only by one additional 
iteration; in the results shown below we employ an new implementation of the methods where more than 
ten additional iterations are included if relevant. A similar idea, where tensor coefficients are iteratively 
determined from higher-rank tensors has been described in Ref. [204] for the massless case. 

(iii) For 5- and 6-point integrals, direct reductions to 5- and 4-point integrals, respectively, are pos- 
sible owing to the four-dimensionality of space-time. For scalar integrals such a reduction was already 
derived in the 1960s [205]. In Refs. [200,202] we follow basically the same strategy to reduce tensor 
integrals, which has the advantage that no inverse Gram determinants appear in the reduction. Instead 
modified Cayley determinants occur in the denominator, but we did not find numerical problems with 
these factors. A reduction similar to ours has been proposed in Ref. [184]. 

We would like to stress two important features of our approach. 

(i) The methods are valid for massive and massless cases. The formulas given in Refs. [200, 202] 
are valid without modifications if IR divergences are regularized with mass parameters or dimension- 
ally. 19 Finite masses can be either real or complex. 

(ii) The in/out structure of the methods is the same as for conventional PV reduction, i.e. no specific 
algebraic manipulations are needed in applications. Therefore, the whole method can be (and in fact is) 
organized as a numerical library for scalar integrals and tensor coefficients. 

We conclude this overview with some comments resulting from our experience collected in the 
treatment of a full 2 — > 4 scattering reaction. 

(i) For a specific point in a multi-particle (multi-parameter) phase space it is highly non-trivial to 
figure out which of the various methods is the most precise. It seems hopeless to split the phase space 
into regions that are dedicated to a given method. Therefore, we estimate the accuracy for the different 
methods at each phase-space point and take the variant promising the highest precision. The accuracy 

"For the method of Ref. [202], this has been shown in Ref. [206]. 
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of the PV method is valued by checking symmetries and PV relations, and by estimating cancellations. 
In the expansion approach, we estimate the number of valid digits based on the expected accuracy of 
the expansions and possible numerical cancellations before the evaluation of the coefficients. In the 
seminumerical approach, the integration error is propagated to the tensor coefficients, together with an 
estimate of possible cancellations. 

(ii) In a complicated phase space it may happen that none of the various methods is perfect or good 
in some exceptional situations. Usually the corresponding events do not significantly contribute to cross 
sections. This issue can only be fathomed in actual applications. To be on the safe side, we employ the 
two independent "rescue systems" with different advantages and limitations. 

(iii) In view of this, figures as shown below are nice illustrations, but should always be taken with 
a grain of salt. No matter how many of such figures are shown, they will never be exhaustive, so that no 
quantitative conclusions on the overall precision of methods can be drawn. 

10.4 Two benchmark phase-space points 

In the following two examples of exceptional phase-space configurations are considered: 20 one with 
small Gram determinant \Z\, another with both \Z\ and \X\ small. These two cases were already qual- 
itatively illustrated in Ref. [207], but without providing explicit numbers. We also note that a complex 
Z-boson mass was used there. Here we switch to a real- valued Z mass to make it easier for other groups to 
compare with our numbers. For the sake of brevity, no results of the seminumerical method are included 
below; such results are illustrated in Ref. [207]. 

10.4.1 A case with a small Gram determinant 

Figure 25 defines the first benchmark point for a 4-point function in which the Gram determinant \Z\ 
becomes small. We compare results of PV reduction with results of the expansion in the small Gram 
determinant as described in Section 5.4 of Ref. [200]. In the upper half of the figure a hexagon diagram 
is shown that contains a box subdiagram with the considered kinematical configuration. The structural 
diagram illustrates the kinematical assignment with internal masses and squared external momenta given 
at the respective lines. The invariants near the arcs are the squares of the sum of momenta flowing into 
the two neighbouring external lines. The explicit values of the masses and invariants are given in the 
figure. As indicated there, the Gram determinant vanishes if the invariant t c ^ approaches the critical 
value i cr it, corresponding to an inner phase-space point. In the plots of Fig. 25 we show results on a 
few tensor coefficients when t c ^ is varied while keeping all other invariants fixed. The variation in i e g 
is translated into a variation of the dimensionless variable x = t e ^/t CI i t — 1 where the exceptional point 
with \Z\ = corresponds to x = 0. 

It is clearly seen in the plot on the l.h.s. that the tensor coefficients calculated with PV reduction 
show numerical instabilities for small x, while the results of the expansion method behave smoothly. The 
PV instabilities increase with increasing tensor rank. The plot on the r.h.s. shows the relative difference 
between the PV results and the corresponding "best" predictions, which are either obtained with the 
PV or the expansion method. With decreasing x this difference rises because of the PV instabilities, 
and for a sufficiently high x the difference becomes zero (and falls out of the plot range), because PV 
reduction promises better accuracy there. It is essential to see a broad region in x where the difference 
is small for each tensor coefficient. This region corresponds to the overlap in which both PV reduction 
and the expansion method are trustworthy, the difference reflecting the uncertainty of the less precise 
result. The plot suggests that both methods should be precise within a relative accuracy of about 10~ 9 
for the considered coefficients which go up to rank 3. As already mentioned for the x values of the 
shown points, the error estimate of the expansion promises better precision, otherwise (for larger x) PV 

20 We have to restrict the set of numerical results to a few selected tensor coefficients; more results can be found under 

http : //wwwth . mppmu . mpg . de /members/ dittmair /tensint s /benchmarks .html . 
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Full diagram: 




Subdiagranr 






Absolute predictions 



Relative deviations from "best" 



fill 



D 



(etc.) Passarino-Vcltman reduction 
— improved by Gram expansion 



10 



10 



10 



10 




Kinematics: 

M z = 91.1876 GeV 
s P ™ = +2xl0 4 GeV 2 
s Pu =+lxl0 4 GeV 2 
fe M =-4xl0 4 GeV 2 

^crit — 

= -6xl0 4 GeV 2 



\Z\^0 & 

x = ^ed/^crit — 1 



Fig. 25: A typical example for 4-point integrals with small \Z\ (x —> 0). The full diagram and the relevant subdiagram are 
given above; absolute predictions (in arbitrary units) for some tensor coefficients, relative deviations from PV reduction, and 
the kinematic specifications are shown below. The precise kinematical assignment is D...(t c ^, sp u , 0, 0, te^, s^u, 0, 0, 0, M|). 
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i 0.18321620694255 
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Table 4: Numerical results corresponding to Fig. 25. 
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Fig. 26: An example for 4-point integrals with both \Z\ and \X\ small (x — + 0). Details as in Fig. 25. The precise kinematical 
assignment is D... (m 2 , s^, 0, s, s mSu , s^d, 0, m 2 , 0, Mf ). 
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1.8810483898149- 
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-1.5970779937221 


fi 67693.244541619 
+ i 0.39078443860164 



Table 5: Numerical results corresponding to Fig. 26. 
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reduction seems to be better. Table 4 provides explicit numbers for the considered tensor coefficients at 
some x values. These numbers could serve as a benchmark also for other methods. 

We recall that the expansion for small \Z\is limited to the case where Xqj and Zjy are not too small 
for at least one set of indices j, k, I. If all Xqj are small, then \X\ is small, too. Such a case is considered 
in the next subsection. The case in which all Z^i are small is elaborated in Section 5.6 of Ref. [200]. 

10.4.2 A case with small Gram and modified Cay ley determinants 

Figure 26 defines the second benchmark point for a 4-point function in which both determinants \Z\ and 
\X\ become small. Here we compare results of PV reduction with results of a simultaneous expansion in 
\Z\ and \X\ as described in Section 5.5 of Ref. [200]. In the upper half of the figure a pentagon diagram 
is shown that contains a box subdiagram with the considered kinematical configuration. The structural 
diagram again illustrates the kinematical situation as in the previous case and the explicit values of the 
masses and invariants are given in the figure. The u-quark mass m u is kept only as regulator of the 
mass singularity, i.e. it is only kept non-zero in the logarithm lnm u , but set to zero otherwise. The 
two determinants \Z\ and \X\ vanish if the two conditions = s and s M p u = s^p are fulfilled. We 
explore the neighbourhood of this exceptional configuration on the specific line parametrized by the 
dimensionless variable x = s^/s — 1 = s^u/s^p — 1, while keeping the internal masses and the 
squares of the external momenta fixed. 

The plot on the l.h.s. again illustrates the instabilities for small x in the PV reduction that become 
more serious for higher tensor ranks, while the results of the expansion method behave smoothly. The 
relative difference between the PV and the corresponding "best" prediction is shown on the r.h.s., re- 
vealing the expected increase for x — > 0. For a sufficiently high x the difference becomes zero, because 
PV reduction is more accurate than the expansion. In the overlap region both methods should be pre- 
cise within a relative accuracy of about 10 -6 for the considered coefficients. Table 5 provides explicit 
numbers for the considered tensor coefficients at some x values. 

The expansion method fails if either all Z^ or all Xij are small. Possible treatments of these 
exceptional cases are also described in Ref. [200]. 
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11. SINGULARITIES IN ONE-LOOP AMPLITUDES FROM THE POINT OF VIEW OF RE- 
DUCTION METHODS 21 

11.1 Introduction 

Obtaining radiative collections requires the evaluation of loop Feynman integrals. The simplest, but also 
the most important, loop integrals are one-loop Feynman integrals. Considerable progress has recently 
been made in developing various approaches for calculating one-loop integrals. Today, at least in prin- 
ciple, it is possible to calculate any of them to arbitrary precision no matter how many external legs 
the corresponding Feynman diagram has. Unfortunately, despite huge development, for a practitioner, 
the calculation of amplitudes up to one-loop contributions is still a difficult task. With the increasing 
complexity of the process under consideration, the number of Feynman diagrams whose contributions 
have to be obtained rises very quickly, as does the complexity of the corresponding one-loop Feynman 
integrals which have to be calculated. Therefore, we are forced to automatize our calculations. Use of 
available automatized algorithms helps tremendously, but the moment when calculations of physically 
relevant processes will demand for practical use unacceptable amounts of computer time and memory 

21 Contributed by: G. Duplancic 
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is not far away. To surpass this problem it is necessary to look for new approaches for calculating one- 
loop amplitudes, but also to implement algorithms in a more "computer friendly" way, which means less 
computer algebra and more numerics. Unfortunately, numerically oriented codes increase our chances 
to face numerical instabilities. This problem is usually connected with the presence of singularities in 
the functions under consideration. It is known that loop Feynamn integrals have rich singularity struc- 
tures. For that reason, it is important to summarize all that is known about the problem as well as to 
share experience from previously performed calculations of amplitudes. Since reduction to the set of 
basic scalar Feynman integrals is at the heart of most methods for calculating Feynman integrals, here 
we discuss singularities from that point of view. Despite different approaches which can be taken, the 
final decomposition of the given Feynman integral, in terms of predefined set of basic integrals, should 
be unique. Therefore, any approach taken to discuss the singularity structure of the final decomposition 
is equally valid. Here the reduction method based on Refs. [151, 153, 208, 209] is used. 

11.2 Definitions and reduction method 

In order to obtain one-loop amplitudes, integrals of the following type are required, 



I N (D-{u l} ) ^ (/ , 2 ) 2 - D/2 



d D l 



lm ' " " I tip 



d U l 



( 2 -) D nf=i [(l + uy-mt+iep 



(108) 
(109) 



The integral i^...„ p (Iq) is a rank P tensor (scalar) one-loop N -point Feynman integral in D-dimensi- 
onal space-time, where are powers of propagators and I + ri (mi) is the momentum (mass) of particle 
propagating along the corresponding internal line. The momentum I is the loop momentum and the r\ 
are linear combinations of external momenta. The scale fj, is the usual dimensional regularization scale 
and the quantity ie (e > 0) represents an infinitesimal imaginary part which ensures causality and, after 
the integration, determines the correct sign of the imaginary part of the logarithms and dilogarithms. It is 
customary to choose the loop momentum in such a way that one of the momenta r.j vanishes. However, 
for general considerations, it is convenient to keep the symmetry of the integral with respect to the indices 
1,..- ,N. 

It can be shown that every tensor one-loop integral can be expressed as a linear combination of 
scalar one-loop integrals by the following equation, 
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where {[(/^[nP 1 ■ ■ ■ [ r N]^ N }m— n P represents a symmetric (with respect to \i\ ■ ■ ■ \ip) combination of 
tensors, each term of which is composed of k metric tensors and ji momenta r,. Therefore, the problem 
of calculating tensor integrals has been reduced to the calculation of the general scalar integral, which is 
the most convenient to evaluate from the following representation, 
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Direct evaluation of the general scalar integral represents a non-trivial problem. However, with the help 
of the recursion relations, the problem can be simplified in the sense that the calculation of the original 
scalar integral can be reduced to the calculation of a certain number of simpler basic integrals. All 
relevant recursion relations for scalar integrals can be written in matrix notation as 
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m 2 . In the following we introduce the notation Sn for the (N + 1) 



where Ri 

(N + 1) matrix in Eq. (112). Making use of relations which follow from Eq.(112), each scalar integral 
Iq (D; {ui}) can be represented as a linear combination of integrals Iq (D'\ {1}) and integrals with the 
number of propagators which is less than N (it has be understood that Iq (D; {■ ■ ■ vi-\, 0, ui + \ • • • }) 
= Iq~ 1 (D; {■ ■ ■ vi-i, vi+i • • • }) )■ For the dimension D' , one usually chooses 4 + 2e, where e is 
the infinitesimal parameter regulating the divergences. By successively applying the above mentioned 
procedure to the remaining less-than-iV-point integrals it is at the end possible to express the integral 
Iq(D; {vi}) as a linear combination of integrals Iq(D'; {1}), k = 1, . . . , N. It is convenient to write 
these basic integrals as 
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where the properties of the 5 function were used. 



11.3 Singularities 

The necessary conditions for Feynman integrals to have singularities are given by the Landau equations. 
In the integral representations given by Eqs. (Ill) and (113), the singularity conditions [192] are given 
by 

JV 

J] I / i (% + 2ie) 1 / i = (114) 



JV 



and 

either y, L = or ^ (Rij + 2ie) pj = for each i. (115) 

Notice that condition (114) is automatically satisfied when conditions (115) are. The singularity of the 
given Feynman integral corresponding to all yi 7^ is called the leading singularity of the integral, while 
those corresponding to some yi = are called lower-order singularities of the integral. Lower-order 
singularities are leading singularities of integrals where all propagators associated with vanishing yjS 
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have been omitted. In the language of Feynman diagrams this translates as contraction to a point of all 
lines associated with = 0. 

Finding the general solution of the Landau equations is non-trivial task. Here we consider only real 
singularities. Real singularities are those occurring for real values of the invariants Rij on the physical 
sheet. Notice that these real values of the invariants do not necessarily correspond to a physically possible 
kinematical configuration. 

Due to presence of the ie, no singularity appears along the real contour of integration in the para- 
metric space in Eqs. (Ill) and (113). It should be understood that singularities appear only in the limit 
ie 0. 

In the previous section, it was described how to express an arbitrary Feynman integral as a linear 
combination of the basic scalar integrals. The question arises if all singularities of the starting integral 
correspond to singularities of the basic scalar integrals or some of them correspond to singularities of 
coefficients of the decomposition. To answer that question, it is enough to check if ie appears in de- 
nominators of the coefficients. That is, singularities appear only in the limit ie —> and if some of 
the coefficients diverge independently of that limit, then the corresponding singularity is artificial in the 
sense that it is not a singularity of the starting Feynman integral. Consequently, such a divergence should 
cancel in sum of all terms in the decomposition. 

The simplest way to see when ie appears in denominators is to invert Eq. (112) by multiplying it 
by inverse of SV- The resulting equation is 
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where 5^ is minor of Sn obtained by removing ith row and jth column, and Rn is an N x N matrix 
with elements equal to Rij. The matrix R^ is sometimes called modified Cayley matrix and its deter- 
minant the modified Cayley determinant. All minors appearing in Eq. (116) as well as Det[SV] are ie 
independent. In all that determinants the first row or column can be simply used to remove completely 
the ie dependance. Therefore, only the recursion relation following from the first row of Eq. (116) will 
have an ie dependent coefficient. From the form of that relation it follows that real singularities can 
appear in the coefficients of decomposition only if relations of that type are used during reduction to 
increase dimension of integrals. In that case, a singularity can appear when Det[i?jy] = 0- As expected, 
the singularity is related to the same matrix which appears in Landau equations (1 14) and (1 15). 

What happens if we calculate the integral exactly for kinematical variables and masses for which 
Det[i?jv] vanishes? In that case the limit ie — ► should produce a divergence. But, from the beginning, 
dimensional regularization was introduced exactly to avoid explicit appearance of divergences. Hence, 
the limit ie — > can be applied and divergences appear in the form of powers of 1/e. It follows that the 
term Det[i?yv] — 2ie Det[5Ar] vanishes and the first row from Eq. (1 16) can be used to reduce the N -point 
integral to a linear combination of (N — 1) -point integrals. 

To complete the discussion, it is necessary to comment on how reduction works for vanishing 
Det[S;v]- Let us first express Det[SV] in a better known form. By subtracting the last column from 
the second, third, . . ., and iVth column, and then the last row from the second, third, . . ., and iVth row, 
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Det[SV] is given by 



Det[S N ] = -Det [-2(r, - r N ) 
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The determinant on the right hand side of Eq. (117) is known as the Gram determinant. If Det[5jv], i.e. 
the Gram determinant, vanishes, then the rows (columns) of the matrix in Eq. (112) should be linearly 
dependent. That is, there are real constants —C, z±, zn, not all of them equal zero, which satisfy the 
equation 
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To see that the constants — C, z\, £jv should be real, just remove the complete ie dependance from 
the system in Eq. (118) by subtracting the equation from the first row multiplied by 2ie from equations 
in all other rows. After multiplying Eq. (1 12) by row ( — C Z\ Z2 ••• zn ), the following relation 
emerges, 
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It is easy to see that by using above relation it is always possible to reduce relevant iV-point scalar integral 
to a linear combination of N — 1-point scalar integrals. For details see [153]. 

From the considerations above, we can conclude that vanishing of the Gram determinant is not 
related to the singularities of Feynman integrals. It is important to point out that the situation is not 
so simple in the case of diagrams with more than one loop. There, Gram determinants are related to 
so-called second-type singularities. 



11.4 Practice and problems 

In practice we deal with 4-dimensional Minkowski space. An immediate consequence of this is that, 
for all integrals with N > 5, Det[SV] vanishes due to the linear dependence of the vectors n and all 
integrals with N > 5 can be reduced to the integrals with N < 5. In view of what has been said above, 
all one-loop integrals are expressible in terms of the integrals 7q (4 + 2e; {1}) with nonvanishing Det[5fc] 
and Det[i?fc], where k = 1, . . . , 5. In fact, for practical calculations, also the 5-point basic scalar integral 
is reducible. That is because we are interested in calculations up to 0{e). Details can be found in the 
literature [150,152,209]. 

For most practical calculations the starting Feynman integrals obtained from Feynman diagrams by 
using Feynman rules are in 4+2e dimensions and with v\ = 1. In the next step, tensor decomposition, Eq. 
(110), will produce scalar integrals with higher dimensions and powers of propagators. By successively 
using all recursion relations following from Eq. (1 16), except the one coming from the first row, in the 
cases of nonvanishing Gram determinants and recursion relations following from Eqs. (118) and (119) 
in the cases of vanishing Gram determinants, it is possible to express an arbitrary Feynman integral 
as a linear combination of integrals Iq(2ti + 2e;{l}) with nonvanishing Det[5fc] and Det[i?^], where 
k = 1, . . . , 5. The possible values for parameter n depend on kinematics involved. If the kinematics 
is such that during reduction no case appears where the constant C in Eqs. (118) and (119) vanishes, 
the parameter n is an integer greater than 1. Now, the recursion from the first row of Eq. (116) can be 
successively used to lower all dimensions down to 4 + 2s. Since in the above procedure that relation 
was never used to increase dimension, from what has been said in the previous section, it follows that 
all singularities are in basic scalar integrals and divergences appearing in coefficients should cancel in 
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the sum. The cases with vanishing C appeal - regularly when dealing with diagrams containing collinear 
external lines, i.e. for exceptional kinematics. 

Assuming the situation described in the previous paragraph, many Gram determinants to different 
powers will appear in denominators of the coefficients when an arbitrary Feynman integral is decomposed 
into the basic integrals. The real problem in practice is when one has to calculate in a kinematical region 
where some of those determinants are small. Since vanishing of the Gram determinant does not corre- 
spond to a singularity, one faces cancellation of big numbers and consequently numerical instabilities. 
In principle, if one is using methods where all Feynman integrals are expressed as linear combinations 
of basic integrals, this problem is unavoidable no matter in which framework coefficients are calculated. 
That is because the decomposition into the basic integrals is unique. However, there are some hints from 
experience as to where one should look to soften this problem. The main guideline is to try to avoid 
separate calculation of diagrams contributing to the process under consideration. Namely, powers of 
determinants in denominators tend to be smaller if a group of diagrams (for example, a gauge invariant 
group) is calculated together. Additionally, one has to use all symmetries of the basic integrals to reduce 
the basic set as much as possible. Of course, at the end, to get more precision, it is always necessary 
to make an expansion around a point where the Gram determinant vanishes. However, if calculating 
in the neighborhood of the point where both Gram and Cayley determinants vanish simultaneously, the 
expansion is problematic because the decomposition is not analytic at that point. One can hope that such 
regions will not give sizable contribution to calculated physical quantities. 

11.5 Conclusion 

Vanishing of various Gram and modified Cayley determinants will always produce numerical instabili- 
ties if reduction methods are used to perform the calculation. The instabilities can be softened by using 
various clever approaches but the question remains, will that work for all practical cases? One can also 
doubt if reduction to basic integrals is the optimal approach to perform calculations which, due to their 
complexity, become more and more numerically oriented. Maybe some kind of direct numerical inte- 
gration of the Feynman integrals is more efficient. Surely this is a more natural approach for numerical 
calculations. 
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Part III 

CROSS SECTIONS 

12. TUNED COMPARISON OF QCD CORRECTIONS TO pp -> WW+jet+X AT THE 
LHC 22 

12.1 Introduction 

The complicated hadron collider environment of the LHC requires not only sufficiently precise predic- 
tions for the expected signals, but also reliable rates for complicated background reactions, especially for 
those that cannot be entirely measured from data. Among such background processes, several involve 
three, four, or even more particles in the final state, rendering the necessary next-to-leading-order (NLO) 
calculations in QCD technically challenging. At the previous Les Houches workshop this problem lead 
to the creation of a list of calculations that are a priority for LHC analyses, the so called "experimenters' 

"Contributed by: T. Binoth, J. Campbell, S. Dittmaier, R.K. Ellis, J.-P. Guillet, S. Kallweit, S. Karg, N. Kauer, G. San- 
guinetti, P. Uwer, G. Zanderighi 
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wishlist for NLO calculations" [16,201]. The process pp — > W + W +jet+X made it to the top of this 
list. 

The process of WW+jet production is an important source for background to the production of a 
Higgs boson that subsequently decays into a W-boson pair, where additional jet activity might arise from 
the production [210]. WW+jet production delivers also potential background to new -physics searches, 
such as the search for supersymmetric particles, because of leptons and missing transverse momentum 
from the W decays. Last, but not least, the process is interesting in its own right, since W-pair production 
processes enable a direct precise analysis of the non-abelian gauge-boson self-interactions, and a large 
fraction of W pairs will show up with additional jet activity at the LHC. 

First results on the calculation of NLO QCD corrections to WW+jet production have been pre- 
sented by two groups in Refs. [3,4]. A third calculation is in progress [211]. In the following the key 
features of these three independent calculations are described and results of an ongoing tuned comparison 
are presented. 

12.2 Descriptions of the various calculations 

At leading order (LO), hadronic WW+jet production receives contributions from the partonic processes 
qq — > W + W~g, qg — ► W + W~q, and qg — ► W + W~q, where q stands for up- or down-type quarks. 
All three channels are related by crossing symmetry. 

The virtual corrections modify the partonic processes that are already present at LO. At NLO 
these corrections are induced by self-energy, vertex, box (4-point), and pentagon (5-point) corrections, 
the latter being the most complicated loop diagrams. Apart from an efficient handling of the huge amount 
of algebra, the most subtle point certainly is the numerically stable evaluation of the numerous tensor 
loop integrals, in particular in the vicinity of exceptional phase-space points. The three calculations 
described below employ completely different loop methods. Some of them are already briefly reviewed 
in Ref. [201], where more details on problems in multi-leg loop calculations and brief descriptions of 
proposed solutions can be found. 

The real corrections are induced by the large variety of processes that result from crossing any 
pair of QCD partons in — > W + W~ qqgg and — ► W + W~ qqq' q' into the initial state. Here the 
main complication in the evaluation is connected to an efficient phase-space integration with a proper 
separation of soft and collinear singularities. For the separation of singularities the three calculations 
all employ the subtraction method [212] using the dipole subtraction formalism of Catani and Seymour 
[175]. 

The calculation of DKU [3] 

This calculation is actually based on two completely independent evaluations of the virtual and real 
corrections. The W bosons are taken to be on shell, but the results on cross sections presented in Ref. [3] 
do not depend on the details of the W decays. 

Both evaluations of loop diagrams start with an amplitude generation by FeynArts, using the two 
independent versions 1.0 [213] and 3.2 [214]. One of the calculations essentially follows the same 
strategy already applied to the related processes of ttH [215] and tt+jet [216] production. Here the 
amplitudes are further processed with in-house Mathematica routines, which automatically create an 
output in Fortran. The IR (soft and collinear) singularities are treated in dimensional regularization and 
analytically separated from the finite remainder as described in Refs. [206,215]. The pentagon tensor 
integrals are directly reduced to box integrals following Ref. [202]. Box and lower-point integrals are 
reduced a la Passarino-Veltman [167] to scalar integrals, which are either calculated analytically or using 
the results of Refs. [148,217,218]. The second loop calculation is based on FormCalc 5.2 [14], which 
automatically produces Fortran code. The reduction of tensor to scalar integrals is done with the help of 
the LoopTools library [14], which also employs the method of Ref. [202] for the 5-point tensor integrals, 
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Passarino-Veltman [167] reduction for the lower-point tensors, and the FF package [219,220] for the 
evaluation of regular scalar integrals. The dimensionally regularized soft or collinear singular 3- and 
4-point integrals had to be added to this library. 

One calculation of the real corrections employs analytical results for helicity amplitudes obtained 
in a spinor formalism. The phase-space integration is performed by a multi-channel Monte Carlo integra- 
tor [221] with weight optimization [222] written in C++. The results for cross sections with two resolved 
hard jets have been checked against results obtained with Whizard 1.50 [223] and Sherpa 1.0.8 [224]. 
Details on this part of the calculation can be found in Ref. [225]. The second evaluation of the real 
corrections is based on scattering amplitudes calculated with Madgraph [158] generated code. The code 
has been modified to allow for a non-diagonal quark mixing matrix and the extraction of the required 
colour and spin structure. The latter enter the evaluation of the dipoles in the Catani-Seymour subtrac- 
tion method. The evaluation of the individual dipoles was performed using a C++ library developed 
during the calculation of the NLO corrections for tt+jet [216]. For the phase-space integration a simple 
mapping has been used where the phase space is generated from a sequential splitting. 

The calculation of CEZ [4] 

The method of choice for calculation of the virtual corrections of Ref. [4] is similar to the techniques 
adopted by the other groups and is based on the semi-numerical method of Ref. [226] augmented with a 
mechanism to handle exceptional configurations [227]. This method has already been used for the NLO 
calculation of Higgs plus dijet production via gluon-gluon fusion [5]. Tree-level matrix elements for real 
radiation have been checked against the results of Madgraph [228]. Soft and collinear singularities are 
handled using the dipole subtraction scheme [175]. As for the other authors, CEZ have performed several 
checks to test the reliability of their code. These include checks of Ward identities of the amplitudes 
containing external gluons. 

The calculation of Ref. [4] is however different from the other two in that the decay of the 
W bosons is included from the outset. Rather than summing over the polarizations of a W boson of 
momentum k with 



the authors of this paper project out the combination of polarizations which occurs in the physical decay 
of the W boson, W~{k) -> e~{h) + v{h), 



The inclusion of the decay is well-motivated from a physical point of view, because it allows phenomeno- 
logical analyses which include cuts on the decay leptons. 

For the purposes of the comparison of virtual matrix elements for a fixed phase-space point, the 
results including the decays can be used to extract the result for the amplitude squared summed over the 
polarization of the vector boson, as would be obtained using Eq. (120). This is achieved by performing 
6 x 6=36 evaluations of the amplitude squared [229] in which each lepton is emitted along three orthog- 
onal axes (in both positive and negative directions) in the corresponding vector-boson center-of-mass 
frame. The results of this comparison, with input parameters tuned for the comparison, will be given 
below. 

The calculation of BGKKS [211] 

This calculation is also done in two independent ways. The graph generation is based on QGRAF [165] 
and was cross checked by having two independent codes. All diagrams neglect the quarks of the 3rd 
generation. 




(120) 




(121) 
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Up to now the LO part and the virtual corrections are evaluated. By using the spinor helicity 
formalism, projectors on the different helicity amplitudes are defined. In this way all Lorentz indices can 
be saturated such that the complexity of the one-loop 5 -point tensor reduction is such that at most rank-1 
5-point integrals appear. For each helicity amplitude an algebraic representation in terms of certain basis 
functions is obtained by using the reduction methods developed in Refs. [152, 184]. The whole algebra 
is done in an automated way by using FORM [166] and MAPLE. In both approaches the IR divergent 
integrals are isolated by using 6-dimensional IR finite box functions such that IR poles are in 3-point 
functions only. One implementation uses the function set defined in Appendix C of Ref. [184], and uses 
the implementation of the Fortran 90 code golem9 0. The other computation uses standard scalar 2- 
and 3-point functions as a basis. The complete algebraic reduction to d=6 scalar box and d=n scalar 2- 
and 3-point functions is largely equivalent to a standard Passarino-Veltman reduction. Only the 5-point 
functions are treated differently [184]. Tractable analytical expressions of the coefficients to the two sets 
of basis functions are obtained for each independent helicity amplitude. 

Discrete symmetries (Bose,C,P) are used to check and relate helicity amplitudes with each other. 
The coefficients are exported to a Fortran code and used to evaluate the loop correction of the process. 

For the treatment of 7,5 the 't Hooft-Veltman scheme is applied. The 7-algebra and the loop 
momenta are split into 4- and (D — 4) -dimensional parts. Whereas the 75 anti-commutes with the D = 4 
matrices, it commutes with the gamma matrices defined in d = D — 4. As is well known the QCD 
corrections of an axial vector current are different from the vector part and a finite renormalisation has 
to be performed. The following gauge boson vertex which includes a finite counterterm for the axial part 
(see e.g. Refs. [230-232]) is used, 

Vy q <i~9v^ + Z b g al ^ with Z 5 = l-C F -±, (122) 

to reinforce the correct chiral structure of the amplitudes. Note that the 't Hooft-Veltman scheme treats 
the observed particles in 4 dimensions but the soft/collinear gluons in D dimensions. This guarantees 
that for the IR subtractions the same Catani-Seymour dipole terms as for conventional dimensional 
regularisation can be used [233]. 

12.3 Tuned comparison of results 

The following results essentially employ the setup of Ref. [3]. The CTEQ6 [28, 234] set of partem 
distribution functions (PDFs) is used throughout, i.e. CTEQ6L1 PDFs with a 1-loop running a s are 
taken in LO and CTEQ6M PDFs with a 2-loop running a s in NLO. Bottom quarks in the initial or final 
states are not included, because the bottom PDF is suppressed w.r.t. to the others. Quark mixing between 
the first two generations is introduced via a Cabibbo angle 9c = 0.227. In the strong coupling constant 
the number of active flavours is Np = 5, and the respective QCD parameters are A^o = 165 MeV and 
Ajp = 226 MeV, leading to o^°{M w ) = 0.13241687663294 and af LO (M w ) = 0.12026290039064. 
The top-quark loop in the gluon self-energy is subtracted at zero momentum. The running of a s is, 
thus, generated solely by the contributions of the light quark and gluon loops. In all results shown in 
the following, the renormalization and factorization scales are set to M\y- The top-quark mass is m t = 
174.3 GeV, the masses of all other quarks are neglected. The weak boson masses are M\y = 80.425 GeV, 
Mz = 91.1876 GeV, and Mh = 150 GeV. The weak mixing angle is set to its on-shell value, i.e. fixed 
by c w = 1 ~~ s w = Myy/Mg, and the electromagnetic coupling constant a is derived from Fermi's 
constant G M = 1.16637 • 10~ 5 GeV~ 2 according to a = VZG^M^sl/ir. 

We apply the jet algorithm of Ref. [235] with R = 1 for the definition of the tagged hard jet and 
restrict the transverse momentum of the hardest jet by prjet > 100 GeV. 
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Table 6: Results for squared LO matrix elements at the phase-space point (123). 

12.4 Results for a single phase-space point 

For the comparison the following set of four-momenta is chosen, 

= (7000,0,0,7000), p£ = (7000,0,0,-7000), (123) 

= (6921.316234371218,3840.577592920205,0,5757.439881432096), 

p% = (772.3825553565997, -67.12960601170266, -279.4421082776151, -712.3990141151700), 

v ^ = (6306.301210272182, -3773.447986908503, 279.4421082776151, -5045.040867316925), 

where the momentum assignment is for a{p\)b(p2) — ► W + (ps)W~ {pi)c{pz) . 

Table 6 shows some results for the (spin- and colour-summed) squared LO matrix elements, as 
obtained with Madgmph [158]. The results of all three groups agree with these numbers within about 13 
digits. 

Because of the different treatment of the number of active flavours in the calculations of DKU and 
CEZ and in order to be independent of the subtraction scheme to cancel IR divergences, we found it useful 
to compare virtual results prior to any subtraction. The 0(a s ) contribution to the virtual, renormalized 
squared amplitude is given by the interference between tree-level and one-loop virtual amplitude, which 
we denote schematically as 

2Re{M£ • M LO } = e\ 2 /(/i rcn ) (c^ + + c ) , (124) 

with 23 /(//ren) = T(l + e ) (4vr/i 2 cn /M^) e and the number of space-time dimensions D = 4 — 2e. In the 
following we split the coefficients of the double and single pole and for the constant part, c_2, c~\, and 
Co, into bosonic contributions ("bos") without closed fermion loops and the remaining fermionic parts. 
The fermionic corrections are further split into contributions from the first two generations ("ferml+2") 
and from the third generation. 

Table 7 shows the results for the bosonic parts of the coefficients c_2, c_i, and Co (c_2 does not 
receive fermion-loop corrections). The results on cq obtained by the different groups typically agree 
within 7—11 digits; the ones on c_2 and c_i agree much better, because they are much easier to calculate. 
The results for the fermionic contributions of the first two generations are given in Table 8. Compared to 
the bosonic corrections these contributions are suppressed by three orders of magnitude. Counting this 
suppression factor, which results from cancellations, as significant digits, the finite parts agree within 
6—9 digits. The agreement is somewhat better in the coefficients of the single pole, which entirely stems 
from the counterterm of the fermion-loop part of the gluon self-energy. The remaining contributions 
from closed loops of the third quark generation are not compared yet. For future reference we show the 
full corrections including all bosonic and fermionic contributions in Table 9. 

23 Note that this factor differs from the overall factor cr extracted when quoting results for one phase-space point in the CEZ 
paper. 
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Table 7: Results for the bosonic virtual corrections at the phase-space point (123) with c-2, c_i and co are defined in Eq. (124). 
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Table 8: Results for the fermionic contributions of the first two quark generations to c_i and Co at the phase-space point (123). 
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Table 9: Results for the full bosonic+fermionic contributions to c„2, c_i and Co at the phase-space point (123). 
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pp^W+W-+jet+X <TL [fb] ciNLo[fb] er^+ilfb] 

DKU 10371.7(12) 14677.6(98) -881.5(42) 

CEZ 10372.26(97) 

BGKKS 10371.7(11) 

Table 10: Results for contributions to the integrated pp cross sections at the LHC in LO and NLO. 

12.5 Results for integrated cross sections 

A tuned comparison of integrated cross sections is still in progress. Table 10 illustrates the agreement 
in the LO cross section obtained by the different groups and provides the DKU result in NLO for future 
comparisons. The subcontribution <7 v i rt +i corresponds to the IR-finite sum of the virtual corrections and 
the contribution of the / operator that is extracted from the real corrections with the dipole subtraction 
formalism [175]. 

12.6 Conclusions 

We have reported on an ongoing tuned comparison of NLO QCD calculations to WW+jet production at 
the LHC. For a fixed phase-space point, the virtual corrections obtained by three different groups using 
different calculational techniques agree within 6-9 digits. The comparison of full NLO cross sections, 
which involve the non-trivial integration of the virtual corrections over the phase space, is still in progress. 

The agreement found so far gives us confidence in the conclusions drawn for physical quantities, 
which were reported in Refs. [3, 4]. 

Acknowledgements 

RU. is supported as Heisenberg Fellow of the Deutsche Forschungsgemeinschaft DFG. This work is 
supported in part by the European Community's Marie-Curie Research Training Network HEPTOOLS 
under contract MRTN-CT-2006-035505 and by the DFG Sonderforschungsbereich/Transregio 9 "Com- 
putergestiitzte Theoretische Teilchenphysik" SFB/TR9. T.B. is supported by the British Science and 
Technology Facilities Council (STFC), the Scottish Universities Physics Alliance (SUPA). T.B., S. Karg 
and N.K. are supported by the Deutsche Forschungsgemeinschaft (DFG, BI- 1050/1). N.K. is supported 
by the Bundesministerium fur Bildung und Forschung (BMBF) under contract 05HT1WWA2. G.Z. is 
supported by the British Science and Technology Facilities Council (STFC). 

13. FROM THE HIGH ENERGY LIMIT OF MASSIVE QCD AMPLITUDES TO THE FULL 
MASS DEPENDENCE 24 

13.1 Introduction 

It is clear that the physics program of the LHC poses new challenges to the theory. In fact, the description 
of hadronic collisions involves several quantities, both non-perturbative and perturbative, the determina- 
tion of which is a highly non-trivial task. As far as the perturbative part is concerned, we are still a 
long way of having the partonic cross sections predicted at a suitable level of accuracy. Whereas most 
processes will have to be known to next-to-leading order, there are some for which the experimental pre- 
cision grants a study going one order higher in the strong coupling constant. Particularly interesting here 
is the top quark pair production cross section. With statistics going into millions of events, a systematics 
dominated error of under 10% is expected already in the first phase of the LHC. Despite years of efforts, 
the appropriate complete NNLO prediction is not yet available. The bottleneck, as in most such cases is 
the evaluation of the two-loop virtual corrections. 

24 Contributed by: M. Czakon 
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Recently, the high energy limit of the amplitudes in the quark annihilation and gluon fusion chan- 
nels has been derived [236, 237] by a mixture of direct evaluation of Feynman graphs and an approach 
based on factorization properties of QCD (see A. Mitov's and S. Moch's contribution). The knowledge 
gained can already be used for the description of high px events and as a test of a future complete pre- 
diction. Clearly due to the behavior of the particle fluxes, what is needed is a calculation covering the 
whole range of variation of the kinematical parameters. It is interesting that one can actually use the high 
energy limit to deal with this problem. Unfortunately, it is not enough to have the whole amplitude, but 
it is rather necessary to know all of the master integrals. In the following, I describe the steps that lead 
to the complete result. 

13.2 The high energy limit 

By the high energy limit, I understand the limit where all the invariants are much larger than the mass. A 
direct approach to the evaluation of the amplitude under this assumption has been devised in [238,239]. 
As a first step, one uses the Laporta algorithm to reduce all of the integrals occurring to a small set of 
masters. In the case at hand, the number of integrals is 145 and 422 for the quark annihilation and gluon 
fusion channels respectively. 

Subsequently, Mellin-Barnes representations are constructed for all the integrals [240,241]. This 
can be done by an automatic package, here by one written by the Author and G. Chachamis 25 . After 
analytic continuation in the dimension of space-time performed with the MB package [243], the integrals 
have the following general form 



where the / function contains, amongst others, a product of T, or possibly ip functions, which have poles 
in z. The desired expansion is obtained by closing the contour and taking residues. As a result, one 
obtains integrals which have lower dimension and a simpler structure. These still require evaluation. 
Due to the fact, that there is a relation between the massive and the massless cases, the result must have a 
similar structure. In particular, it has to be given by harmonic poly logarithms, and therefore it should be 
possible to resum the integrals by further closing contours and evaluating the resulting series. This can 
again be achieved automatically with the help of the XSummer package [244]. 

What remains at the end are integrals, which are pure numbers, but do not have a structure sug- 
gesting a solution in terms of harmonic series. The same argument as before shows, however, that this 
must be the case. Instead of working out specific methods for particular integrals, it turned out to be 
possible to evaluate them to very high precision and subsequently use the PSLQ algorithm to reconstruct 
the solution in terms of Riemann zeta values. 

It has to be noted, that the procedure sketched above works for the majority of cases, but some 
remain at the end. For these, it is usually necessary to change the basis of integrals, in order to obtain 
expressions of suitable structure and/or size for evaluation. At present this program has been completed 
for the quark annihilation channel, and thus all the color structures given in [236] have been computed 
directly with agreement with the factorization approach. The gluon fusion channel is still under way. 

13.3 Power corrections 

As explained in the introduction, the high energy limit by itself is not enough for practical applications. 
To go one step further, it is possible to compute power corrections in the mass. These will then cover 
most of the range, apart from the threshold region and the small angle region, where the series is not 
convergent any more. 

25 There is a public package available [242] that constructs representations for planar graphs. In the present case, also non- 
planar graphs occur, however. 
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Fig. 27: Bare leading color amplitude for top quark pair production in the quark annihilation channel expanded in the mass. 
The more divergent terms at threshold (left of the plots) correspond to higher orders of expansion. The left panel corresponds 
to 90 degree scattering, whereas the right to forward scattering. The variables are defined in the text. 



The main idea is as follows. The derivative of any Feynman integral with respect to any kinemat- 
ical variable is again a Feynman integral with possibly higher powers of denominators or numerators. 
These can, however, be reduced to the same master integrals. This means that one can construct a par- 
tially triangular system of differential equations in the mass [245,246], which can subsequently be solved 
in the form of a power series. 

In Fig. 27, 1 show the result of expansion for the leading color term. The kinematic variable x is 

x = —, t = (p 3 -pi) 2 - m 2 , (126) 



and its variation within the range [1/2(1 — /?), 1/2(1 + (3)], where (3 = yl — 4m| / s is the velocity, 
corresponds to angular variation between the forward and backward scattering. 

The series appears to be asymptotic at the boundaries. Unfortunately, the behavior is worse for the 
subleading color terms, as a consequence of the Coulomb singularity among others. 



13.4 Numerical evaluation 

Using the same system of differential equations one can obtain a full numerical solution to the problem. 
The only requirement is to have the boundary conditions to suitable accuracy. These are provided by 
the series expansions of the previous section. It is crucial to perform the numerical integration along a 
contour in the complex plane, since there are spurious singularities along the real axis. Here, I chose 
an ellipse, because of the improved control on the integration error that one gets from the software used 
(ODEPACK). 

Fig. 28 shows the solution in the range, where the expansion of the previous section starts to di- 
verge. The achievable precision, if double precision arithmetic is used, is about 10 digits for most points, 
with evaluation times of the order of a second. This is going to be substantially slower, when subleading 
color terms will be added. However, the method is fast and precise enough to be sufficient for practi- 
cal applications. In particular it is possible to construct grids of solutions, which will be subsequently 
interpolated when implemented as part of a Monte Carlo program. 

It is clear that the method is suitable for problems, which have a relatively small number of scales, 
and seems to be perfect for 2 — ► 2 QCD processes at the two-loop level. The main drawback is the size 
of the expressions, and the difficulties connected to the derivation of the boundary conditions. 
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13.5 Conclusions 

I have described an approach for the evaluation of massive QCD amplitudes starting from the high energy 
limit and its application to the NNLO corrections to the top quark pair production cross section. Needless 
to say, the same procedure can be applied to other problems of interest. At present the Author, together 
with G. Chachamis and D. Eiras, is working on the corrections to gauge boson pair production. 
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14. MUCH CAN BE SAID ABOUT MASSIVE AMPLITUDES JUST FROM KNOWING 
THEIR MASSLESS LIMIT 26 

14.1 The high-energy limit 

For the precise evaluation of collider observables the knowledge of the pure virtual correction to the 
corresponding Born process is required. This is true at any order (NLO, NNLO, etc). In presence of 
heavy flavors, and especially at higher orders, the problem of their evaluation becomes acute. 

There are very important applications awaiting such results. An example of central importance 
is top production at LHC which is one of the few eagerly awaited precision observables at this collider. 
One of the peculiar features of top production at LHC, and in contrast to the situation at the Tevatron, is 
that no specific kinematical region dominates the cross-section. This is due to the shape of the luminosity 
function for LHC kinematics. 

Direct calculation of the amplitudes is certainly a very demanding task and it seems that one can 
hope that numerical results in some, hopefully easy to handle form, will become available soon (see 
the contribution by M. Czakon for progress in this direction). Here we consider an alternative approach 
which explores the special properties of the gauge theory amplitudes in the high-energy limit and easily 
provides (partial) results for the heavy flavor amplitudes even at higher perturbative orders. 

26 Contributed by: A. Mitov, S. Moch 
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In the following we start by introducing the concept of high-energy limit with the help of simple 
and physically motivated arguments. By high-energy limit one means a kinematical situation where 
the corresponding invariants are much larger than the masses of the heavy particles of interest. In the 
following we will consider the case of a single massive fermion with mass m in presence of a typical 
large kinematical invariant Q. Specific examples are detailed in section 14.3. If the quantity of interest 
(like total or differential cross-section, amplitude, etc.) is regular in the limit m/Q — > then the high- 
energy limit is quite trivial: it is an m-independent function of the kinematical invariants which coincides 
with the one evaluated in the massless limit. Therefore it can be computed by setting the mass m to zero 
from the very beginning. 

Such a situation is, however, relatively rare. In most quantities of interest, like the differential 
ones, the limit m /Q — ► is singular. The obvious manifestation of that singularity in the results is the 
presence of terms of the type ~ In™ (m/Q). When such contributions appear (and in fact this is the 
typical situation) the high-energy limit is defined as the full result with all power corrections in the mass 
neglected, i.e. it contains all logarithms (not multiplied by powers of the mass) as well as the so-called 
"constant" or mass-independent terms. Clearly, in such cases the high-energy limit is different from the 
massless limit. 

Before we detail the relation between these two limits, we would first like to clarify the origin 
and meaning of the logarithmic terms mentioned above. These terms are known as (quasi-) collinear 
logs since they originate from emissions of collinear radiation. To be precise, the role of the mass is to 
regulate small angle emissions that would otherwise diverge in a massless theory; see Ref. [247] for a 
detailed exposition. In this regard, a parton's mass gets dual significance, since one can take one of the 
following two viewpoints: 

• small or large, the mass is nevertheless non-zero, therefore the result is always (collinearly) 

finite; 

• the small mass is just a formal regulator for collinear singularities much like dimensional regu- 
larization in the purely massless case. 

In this write-up we take the unifying viewpoint that both approaches are useful and do not have 
to be considered as alternatives to each other. One can think of the mass as a regulator which is helpful 
in deriving certain properties of the theory but it can also be thought of as an approximation to the full 
massive result which is surprisingly good in many physical applications. 

The prominent role these logarithmic terms play in physical applications has been acknowledged 
long ago, and has been formalized in the so-called Perturbative Fragmentation Function approach [248] 
now known through two-loops [249, 250]; for a recent review see [25 1], The idea behind this formalism 
is the fact that up to power collections in the mass, a differential with respect to some kinematical 
parameter z cross-section for the production of a massive partem h, can be written as: 

^(z, Q,m)=J2 ^r(z, Q) ® D a ^ h {z, m) + 0(m) . (127) 

a 

The function D a _^(z, m) does not depend on the hard scale Q and is thus a process independent 
object that can be computed to any fixed order. It has the important property that it contains all the mass 
dependence within the approximation indicated in Eq. (127). On the other side the partonic cross-section 
do a for the production of any parton a is intrinsically massless, i.e. it is obtained from a calculation where 
m = is set from the very beginning. Of course, collinear singularities are still present in a massless 
calculation but they are regulated dimensionally, i.e. they appeal - as poles in e, where d = 4 — 2e: 

^(z, Q,e)=Y^ ^(z, Q) ® T ba (z, e) . (128) 

b 

The explicit expression for the collinear counterterm T contains arbitrariness; the only condition on it 
is that it contains all poles in e. It has become a standard practice in recent years to work in the MS 
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scheme where T contains only poles. The choice of a subtraction scheme is of course also implicit in the 
definition of the function D in Eq. (127). 

From Eqs. (127), (128) it is quite clear that one can obtain a massive cross-section in the small- 
mass (or high-energy) limit by performing a purely massless calculation. The calculational simplifi- 
cations following from this can be enormous, especially at higher orders. The usefulness of such an 
approach has been appreciated in the past in many applications related to heavy quark production (typ- 
ically b and c) at special kinematics like large Pt hadroproduction and e + e~ annihilation at the Z-pole 
(see [252] for a review). In such kinematical configurations the neglected power corrections can be as 
low as a few percent effect and are often totally negligible. 

A second virtue of Eq. (127), and one that cannot be matched in conventional perturbation theory, 
is that it allows resummation of large collinear logs hx(Q/m) to all perturbative orders. This feature 
is due to the fact that the function D satisfies the DGLAP evolution equation, or in other words one 
achieves exponentiation of the (remnants of) soft and collinear singularities. 

As we will demonstrate in the next section, all these features of massive cross-sections in the 
small-mass limit can be translated to massive amplitudes in gauge theories where similar properties can 
be uncovered. Moreover, one can exploit these properties in much the same way; this is illustrated by the 
physical applications we consider in section 14.3. 

14.2 Factorization in massive amplitudes 

As was indicated above, in the following we will be concerned with the factorization properties of mas- 
sive QCD amplitudes in the high-energy limit. Since one of our main objectives is to relate the small- 
mass limit of an amplitude with its massless limit, we start our discussion with a brief review of the 
well-known factorization properties of massless amplitudes [253,254]. 

The scattering amplitude A4 P 

\M P ) = M p {{h},{m i },{c l }^,a s {ii 2 )^ , (129) 

for a general 2 — > n scattering processes of on-shell partons pi 

p: P1+P2 -> P3~\ VVn+2- (130) 

with a set of fixed external momenta {h}, masses {mj} and color quantum numbers {q}, can be written 
in the massless case m{ = as a product of three functions Jp~°\ 5p m_0 ^ and 

\ Mp) (m=0) =J (m=0) ^^2)^^=0) f {k . h & tas ^2 ) e \ | Wp)? (131) 

The decomposition Eq. (131) can be understood with simple physical arguments. The jet func- 
tion J p m ~^ contains all collinearly sensitive contributions, is color-diagonal and depends only on the 
external partons. On the other side the soft function S p contains all soft radiation interferences and 
is therefore process specific. Finally, the short-distance dynamics of the hard scattering is described by 
the (infrared finite) hard function H p . To leading order this function is just proportional to the Born 
amplitude. More details about the above expressions can be found in the review [255]. 

As was explained in [254], the decomposition Eq. (131) contains arbitrariness related to subleading 
soft as well as finite contributions, which can be removed by fixing a prescription. A convenient and 
natural choice is to identify the jet function with the massless form factor for the flavor corresponding to 
any particular leg, i.e.: 

j(m= 0)= yi 4 m=o) = n (4r 0) )"' (i32) 

i£ {all legs} ig {all legs} 
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where i = q,g for quarks and gluons. J'X ~ is the individual jet function of each external parton. The 
needed massless jet factors are known through three-loops and the soft functions through two-loops for 
any 2 — > n scattering process [256]. 

We are now ready to consider the massive case. Based on our discussion in the previous sections 
in the small-mass limit one should expect a decomposition of massive amplitudes similar to the one in 
Eq. (131). Let us be more specific: we know that in the massive case collinear logs do appear but we 
also know that they should be absorbed in a corresponding jet function. On the other side, up to power 
corrections, the soft and hard functions in the massive case should be the same as in the massless case 
since, by construction, they are not sensitive to collinear emissions. With the exception of contributions 
related to heavy quark loops (to be discussed below) in the presence of a hard scale Q we write for the 
massive amplitudes (130): 

lMp) (m) =J (m) ^, {m .} )aB ( M 2) j ^ 5 (m=o) ({ ki }, a s (^ 2 ), e) \H p ) + 0(m) . (133) 

It is very easy to find out what the jet function in the massive case should be. Working in the 
prescription chosen for the massless case, one can apply the amplitude decomposition to the form factor 
itself; the latter has no nontrivial soft or hard functions. Therefore, in the massive case the jet function 
must be nothing but the massive form factor evaluated in the small-mass limit. 

Combining Eqs. (131), (133) one gets the following very suggestive relation [257]: 

H m) = II ( Z ^f X M^=V+0(m,...), (134) 

i£ {all legs} 

where, 

7 {m\o)(m 2 \ ( m ) / Q 2 m 2 \f ( m =o) f Q 2 W 1 ^ n<xv\ 

z w \^fl ,aB,e ) = k \Jj^ , ^jF ,as,e ) \V\ VT^'VJ (135) 

is a universal, process independent factor. It is sensitive to the definition of the mass m as well as the 
coupling constant (see [257] for details). The process-independence in Eq. (135) is manifest because 
Z^T^ is only a function of the process-independent ratio of scales fi 2 /m 2 . The process-dependent scale 
Q cancels completely between the massive and the massless form factors. 

The last statement, however, requires one important clarification. From the explicit results for the 
massive and massless form factors one can easily see that starting from two loops the ratio indicated 
above contains also Q-dependent logarithmic terms originating from diagrams with the heavy parton in 
loops. It is these terms that we have indicated with dots in Eq. (135). Luckily, these terms are easy to 
recognize and to separate since in the color decomposition of the amplitudes they are proportional to the 
number of heavy flavors n^. For that reason in the definition of Z-factor given originally in Ref. [257] 
contributions proportional to the number of heavy flavors have been excluded, as indicated by the dots 
in Eq. (134). A first step in the understanding of the loop contributions and their incorporation into 
the factorization approach was made in Ref. [258] in the context of Bhabha scattering. We discuss this 
process as well as other applications in the next section. 

Comparing the results of this section with the ones in the previous section, we can clearly see the 
similarities offered by QCD factorization between small-mass limits of amplitudes and cross-sections. 
In both cases the small-mass results are proportional to the corresponding massless results. The pro- 
portionality factors are process independent and universal. The proportionality is in the sense of usual 
multiplication for amplitudes and convolution for cross-sections, as usual. Moreover, it was explained in 
Ref. [257] the so-called Z-factor in Eq. (135) seems in fact to equal the pure virtual contributions to the 
perturbative fragmentation function D in Eq. (127). 
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14.3 Applications 

The results in the previous section have been cross-checked at the amplitude level with the general-mass 
predictions for the structure of the e-poles and ln(m) terms of any one-loop amplitude [247]. Complete 
agreement was found. We have also checked that for the process qq — > hh the prediction based on 
Eq. (134) completely agrees with the results from the one-loop calculation of Ref. [259]. We want to 
stress that we have compared not only the singular terms but also all terms that are finite in the limits 
e — ► and m — > 0. The agreement applies to all color structures of the amplitude as well as for both its 
real and imaginary parts. 

In subsequent work [236, 237] a prediction for the small-mass limit of all two-loop heavy quark 
production squared amplitudes at hadron colliders has been made, while the terms proportional to n h 
were obtained from a direct calculation. We will not go into details here (they can be found for example 
in the recent review [255]) but will only summarize the main features of the result: several of the color 
structures were calculated both directly as well as predicted and we observed full agreement between the 
two approaches. Therefore, this is a first two-loop check for the factorization approach and represents a 
direct confirmation of its validity. 

Another obvious application where the small-mass limit plays important role is Bhabha scattering. 
The knowledge of the two-loop QED massive amplitudes in the small-mass limit there is needed for 
achieving the intended precision of the luminosity measurement; see for example [260]. Complete results 
for the photonic corrections to large-angle Bhabha scattering were first obtained by Penin [261] and later 
confirmed in Ref. [258] in the approach discussed in the previous section. Therefore, this is yet another 
example of its usefulness and power. 

14.4 Conclusions 

We have presented a newly developed relation between massive and massless QCD amplitudes. We 
have emphasized its relevance for physical applications and its ability to seamlessly produce results for 
processes that cannot be calculated currently by direct means. 

The relation was introduced based on the idea for massless limit of a massive amplitude and was 
given in parallel to the much better known relation between massive and massless differential cross- 
sections. 

The new relation between massive and massless amplitudes represents the proper generalization 
of the naive textbook replacement relation 1/e — ► ln(m) + . . . to all perturbative orders and for any 
process. Moreover, with the obvious identification of the color factors, the relation is applicable to any 
SU(N) gauge theory, QCD being a prominent example. QCD and QED applications like heavy quark 
production at hadron colliders at two loops and two-loop corrections to Bhabha scattering were briefly 
discussed. 

15. NNLO PREDICTIONS FOR HADRONIC EVENT SHAPES IN e+e~ ANNIHILATIONS 27 
15.1 Introduction 

For more than a decade experiments at LEP (CERN) and SLC (SLAC) gathered a wealth of high pre- 
cision high energy hadronic data from electron-positron annihilation at a range of centre-of-mass ener- 
gies [262-278]. This data provides one of the cleanest ways of probing our quantitative understanding of 
QCD. This is particularly so because the strong interactions occur only in the final state and are not en- 
tangled with the parton density functions associated with beams of hadrons. As the understanding of the 
strong interaction, and the capability of making more precise theoretical predictions, develops, more and 
more stringent comparisons of theory and experiment are possible, leading to improved measurements 
of fundamental quantities such as the strong coupling constant [252,279]. 

27 Contributed by: G. Dissertori, A. Gehrmann-De Ridder, T. Gehrmann, E.W.N. Glover, G. Heinrich, H. Stenzel 
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In addition to measuring multi-jet production rates, more specific information about the topol- 
ogy of the events can be extracted. To this end, many variables have been introduced which charac- 
terise the hadronic structure of an event. With the precision data from LEP and SLC, experimental 
distributions for such event shape variables have been extensively studied and have been compared 
with theoretical calculations based on next-to-leading order (NLO) parton-level event generator pro- 
grams [212,280-285], improved by resumming kinematically-dominant leading and next-to-leading log- 
arithms (NLO+NLL) [286-291] and by the inclusion of non-perturbative models of power-suppressed 
hadronisation effects [292-295]. 

Comparing the different sources of error in the extraction of a s from hadronic data, one finds that 
the purely experimental error is negligible compared to the theoretical uncertainty. There are two sources 
of theoretical uncertainty: the theoretical description of the parton-to-hadron transition (hadronisation 
uncertainty) and the uncertainty stemming from the truncation of the perturbative series at a certain 
order, as estimated by scale variations (perturbative or scale uncertainty). Although the precise size 
of the hadronisation uncertainty is debatable and perhaps often underestimated, it is conventional to 
consider the scale uncertainty as the dominant source of theoretical error on the precise determination 
of a s from three-jet observables. This scale uncertainty can be lowered only by including perturbative 
QCD corrections beyond NLO. 

We report here on the computation of NNLO corrections to event shape distributions, and discuss 
the impact of these corrections on the extraction of a s from LEP data. 



15.2 Event shape variables 

In order to characterise hadronic final states in electron-positron annihilation, a variety of event shape 
variables have been proposed in the literature, for a review see e.g. [291, 296]. These variables can be 
categorised into different classes, according to the minimal number of final-state particles required for 
them to be non- vanishing: In the following we shall only consider three particle final states which are 
thus closely related to three -jet final states. 

Among those shape variables, six variables were studied in great detail: the thrust T [297,298], 
the normalised heavy jet mass p [299], the wide and total jet broadenings B\v and Bt [300], the C- 
parameter [301, 302] and the transition from three-jet to two-jet final states in the Durham jet algorithm 
y 3 [303-307]. 

The perturbative expansion for the distribution of a generic observable y up to NNLO at e + e~ 
centre-of-mass energy y/s, for a renormalisation scale p 2 , is given by 

1 da, 2 . /a s (fi 2 )\dA faJfi 2 )\ 2 fdB di , fi 2 





cr had dy \ 2vr J dy \ 2vr J \dy dy s 

dC dB n , 
— + 2— A)log^ 
dy dy s 

+^ (/5 2 log 2 - + A log ^) | + G(at) . (136) 

The dimensionless perturbative coefficients A, B and C depend only on the event shape variable y. They 
are computed by a fixed-order parton-level calculation, which includes final states with three partons at 
LO, up to four partons at NLO and up to five partons at NNLO. LO and NLO corrections to event shapes 
have been available already for a long time [212,280-285]. 

The calculation of the NNLO corrections is carried out using a newly developed parton-level event 
generator programme EE RAD 3 which contains the relevant matrix elements with up to five external 
partons [139, 160,308-315]. Besides explicit infrared divergences from the loop integrals, the four- 
parton and five-parton contributions yield infrared divergent contributions if one or two of the final state 
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partons become collinear or soft. In order to extract these infrared divergences and combine them with the 
virtual corrections, the antenna subtraction method [316-318] was extended to NNLO level [319-322] 
and implemented for e + e~ — ► 3 jets and related event-shape variables [323]. The analytical cancellation 
of all infrared divergences serves as a very strong check on the implementation. EERAD3 yields the 
perturbative A, B and C coefficients 28 as histograms for all infrared-safe event-shape variables related 
to three-particle final states at leading order. As a cross check, the A and B coefficients have also 
been obtained from an independent integration [283-285] of the NLO matrix elements [212], showing 
excellent agreement. 

For small values of the event shape variable y, the fixed-order expansion, eq. (136), fails to con- 
verge, because the fixed-order coefficients are enhanced by powers of ln(l/ y). In order to obtain reliable 
predictions in the region of y -C 1 it is necessary to resum entire sets of logarithmic terms at all orders in 
a s . A detailed description of the predictions at next-to-leading-logarithmic approximation (NLLA) can 
be found in Ref. [325]. 

15.3 Generic features of the NNLO corrections 

The precise size and shape of the NNLO corrections depend on the observable in question. Common to 
all observables is the divergent behaviour of the fixed-order prediction in the two-jet limit, where soft- 
gluon effects at all orders become important, and where resummation is needed. For several event shape 
variables (especially T and C) the full kinematical range is not yet covered for three partons, but attained 
only in the multi-jet limit. In this case, the fixed-order description is also not applicable since it is limited 
to a fixed multiplicity (five partons at NNLO). Consequently, the fixed-order description is expected to 
be reliable in a restricted interval bounded by the two-jet limit on one side and the multi-jet limit on the 
other side. 

In this intermediate region, we observe that inclusion of NNLO corrections (evaluated at the Z- 
boson mass, and for a fixed value of the strong coupling constant) typically increases the previously 
available NLO prediction. The magnitude of this increase differs considerably between different observ- 
ables [324, 326], it is substantial for T (18%), B T (17%) and C (15%), moderate for p and B w (both 
10%) and small for I3 (6%). For all shape variables, we observe that the renormalisation scale uncer- 
tainty of the NNLO prediction is reduced by a factor of two or more compared to the NLO prediction. 
Inclusion of the NNLO corrections also modifies the shape of the event shape distributions. We observe 
that the NNLO prediction describes the shape of the measured event shape distributions over a wider 
kinematical range than the NLO prediction, both towards the two-jet and the multi-jet limit. To illus- 
trate the impact of the NNLO corrections, we compare the fixed-order predictions for I3 to LEP2-data 
obtained by the ALPEH experiment in Figure 29, which illustrates especially the improvement when 
approaching the two-jet region, corresponding to large —\n{Y^). 

15.4 Determination of the strong coupling constant 

Event shape data from LEP and LEP2 were used in the past for a precise determination of the strong 
coupling constant a s . These studies were based on the previously available NLO results, improved by 
NLLA resummation; the resulting error on a s was completely dominated by the renormalisation scale 
uncertainty inherent to the NLO calculation. Using the newly computed NNLO corrections to event 
shape variables, we performed a new extraction of a s from data on the standard set of six event shape 
variables, measured by the ALEPH collaboration [263] at centre-of-mass energies of 91.2, 133, 161, 
172, 183, 189, 200 and 206 GeV. The event-shape distributions were obtained using the reconstructed 
momenta and energies of charged and neutral particles. The measurements have been corrected for 
detector effects, ie., the final distributions correspond to a so-called particle (or hadron) level (stable 
hadrons and leptons after hadronisation). 

28 A, B and C differ from A, B and C in their normalisation to ctq instead of tr^ad [324]. 
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-ln(Y 3 ) 

Fig. 29: Perturbative fixed-order predictions for the Y3 -distribution. 

The coupling constant a s is determined from a fit of the perturbative QCD predictions to measured 
event-shape distributions. The procedure adopted here follows closely the one described in Ref. [263]. 
Event-shape distributions are fitted in a central region of the three-jet production, where a good pertur- 
bative description is available. The fit range is placed inside the region where hadronisation and detector 
corrections are below 25% and the signal-to-background ratio at LEP2 is above one. At the higher LEP2 
energies the good perturbative description extends further into the two-jet region, while in the four-jet 
region the background becomes large. Thus the fit range is selected as a result of an iterative procedure 
balancing theoretical, experimental and statistical uncertainties. 

Here we concentrate on fits of NNLO predictions [327] and compare them to pure NLO and 
matched NLO+NLLA predictions as used in the analysis of Ref. [263]. Results from individual event 
shapes are displayed in Figure 30. The combination of all NNLO determinations from all shape variables 
yields 

a,(M§) = 0.1240 ± 0.0008 (stat) ± 0.0010 (exp) ± 0.0011 (had) ± 0.0029 (theo), (137) 

which is indicated by the error band on Figure 30. We observe a clear improvement in the fit quality 
when going to NNLO accuracy. Compared to NLO the value of a s is lowered by about 10%, but still 
higher than for NLO+NLLA [263], which shows the obvious need for a matching of NNLO+NLLA 
for a fully reliable result. The scatter among the a s -values extracted from different shape variables is 
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Fig. 30: The measurements of the strong coupling constant a s for the six event shapes, at s/s = Mz, when using 
QCD predictions at different approximations in perturbation theory. 



lowered considerably, and the theoretical uncertainty is decreased by a factor 2 (1.3) compared to NLO 
(NLO+NLLA). 

These observations visibly illustrate the improvements gained from the inclusion of the NNLO 
corrections, and highlight the need for further studies on the matching of NNLO+NLLA, and on the 
derivation of NNLLA resummation terms. 
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Part IV 

PARTON SHOWERS 

16. DEVELOPMENTS IN LEADING ORDER PARTON SHOWERS 29 

At the Les Houches workshop, there was lively discussion of parton showers as represented in Monte 
Carlo event generators. One of the main current issues in this field is the problem of matrix-element 

29 Contributed by: D.E. Soper, P.Z. Skands 
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/ parton-shower matching, and the workshop saw several reports on ways of (re-)formulating partem 
showers that could make this problem easier to deal with, a trend one might denote with a fancy word 
as "designer showers". In this section, we review the part of that discussion that relates to how a leading 
order parton shower can be organized. Despite the apparent differences, all the new approaches can be 
discussed at a common footing if we adopt a little bit of notation (adapted from [328]). 

A typical parton shower algorithm for hadron-hadron collisions works with states with two ini- 
tial state partons, a and b, and some number m of final state partons that we can label with integers 
1,2, ... ,m. The momenta of these partons can then be specified by giving {p} rn = {pa,,Pb,Pi , ■ ■ ■ , Pm}- 
Each parton also carries a flavor / £ {g, u, u, d, d, . . . }, so that the momenta and flavors can be speci- 
fied with {p, f} rn . Typically, we also keep track of color connections (the labels of one or two partons 
to which parton i is connected in the leading-color limit). We may therefore denote the complete set of 
m + 2 partons by {p, f, c} m , where c denotes the color connections. 

We can now consider the states \{p, f, c} m ) to form a basis for a vector space in the sense of 
statistical mechanics. After some amount of shower evolution starting from a basis state \{p, /, 0)2) 
with two final state partons, one reaches a state |p) that is a linear combination of basis states, so that 
({p, f, c} m \p) represents the probability, in the shower model, for the state \p) to consist of m + 2 
partons with momenta, flavors, and colors {p, f, c} rn . 

As the state develops, partons split. The evolution of the state is tracked with a shower "time" 
t which can be interpreted as (the logarithm of) a typical time for a quantum process such as a parton 
splitting. In most of the current algorithms, the shower time is the logarithm of the virtuality or transverse 
momentum in a splitting. (In Herwig, the shower time represents the energy of the mother parton times 
the square of the splitting angle. In order to cast HERWIG into the form presented here, one also needs 
a cut on virtuality such that splittings with too small virtuality are not allowed. In other parton shower 
algorithms, there is also a smallest virtuality allowed, but that can be obtained by simply stopping the 
shower evolution at some point.) 

The evolution starts with the hard process and works forward in physical time for final state evo- 
lution and backwards in physical time for the evolution of the initial state. Thus we take the shower time 
for a splitting {p, f, c} m — ► {p, f, c} m+ i to be t = t({p, f, c} m+ i) where, for instance if i and j are the 
daughter partons and Q represents the virtuality scale for the hard process that starts the shower, 

K{P, f, g}m+i) = log ( ^ ^ )2 ) . (138) 

It can, and should, be debated whether there is a preferred choice for the shower evolution variable and, 
if so, what it is. 

Using this notation, we can represent what a typical parton shower Monte Carlo does. This repre- 
sentation is an approximation to what real computer codes do. We assume that each stage of evolution 
is independent of what happened at previous stages, depending instead only on the shower time t and 
the partonic state at that stage. This is not the case if, for instance, we do not exactly conserve four- 
momentum at each stage and then adjust the parton four-momenta at the end. 

If we start with a particular basis state \{p, f, c} m ) at shower time to, then at a later time if we get 
a state related to \{p, f, c} m ) by an evolution operator U(t' , to). In the notation of conventional parton 
showers, based on collinear DGLAP splitting kernels, the form of the evolution operator would be 

U {t 1 , t ) j {p, f, c} m ) = A(t', t ; {p, f, c} m )\ {p, f, c} m ) 

+ Y, tdhA^to-ApJ^U) f^di f^^P i ^ jk (z)U(t',t 1 )\{pJ,c} m+1 ) , (139) 

where dt\ = dQ 2 /Q 2 is the differential of the evolution variable, z is an energy-momentum sharing 
fraction, P(z) are the DGLAP splitting kernels, and we include an integral over angle that is usually uni- 
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formly distributed. Once the algorithm picks which partem splits, the flavors, and the splitting variables 
t, z, cj>, the new state {p, f, c} m+ i is known. Reformulating Eq. (139) in the notation outlined above, the 
second term changes appearance slightly, 



U{t',t Q )\{p, f,c} m ) = A(t',t ;{p,f,c} rn )\{p,f,c} m ) 
+ / dh A(ti,t ;{p,f,c} m ) 

Jtn 



x / ({p,/,g}m+i|^l(ti)|{p,/,c} TO ) U{t'M) \{pj,c} m+1 ) . 

(140) 

In either notation, the second term represents that at a shower time t\ > to, the first splitting occurs. 
This splitting time is determined on a probabilistic basis, so t\ is integrated over. The probability to get 
a particular state {p, f, c} m+ i is given by 

({p,/,fiWi|Hi(ti)|{p,/,c} OT ) , (141) 

where Hi is the splitting operator, analogous to the interaction hamiltonian in quantum mechanics. There 
is an integration over the possible outcomes {p, f, c} m+ i. The requirement that the splitting at shower 
time t\ be the first after to means that we must include the probability that there is no earlier splitting. 
This "no-branching" probability is given by a function (the Sudakov form factor) 

A(tutoi{pJ,c} m ) . (142) 

In a lowest order shower, this function is fixed so that the probability not to split in shower time interval 
dt\ is 1 minus the probability to split, 

A(ti,t ;{p,f,c} m ) = 

( f* 1 1 f - \\ ( 143 ) 

exp^-y dr ^— j- ^ J [d{p,f,c} rn+1 \ {{p, /, c} m+ x\Hi{T)\{p, f, c} m ) J . 

The last ingredient in line two of Eq. (140) is the evolution operator U(t', ti). This says that further 
splittings can happen, in the same way, once the first splitting has occurred. It can also happen that there 
is no splitting generated between shower times to and t\. This is represented in the first term of Eq. (140). 

Evidently, the main content of a partem shower resides in the generator TL\(t) of the evolution. 
This has two main parts: a splitting function and a momentum mapping. 

Consider first the splitting functions, functions of the daughter partem momenta that give the prob- 
ability to split. If a parton splits into two nearly collinear partons, then the splitting function must match 
the probability given by Feynman graphs in the collinear limit. For the moment, we discuss a spin 
averaged, leading color shower. Then the splitting function matches the result from Feynman graphs 
averaged over the mother parton spin and summed over the daughter spins in the approximation of ne- 
glecting contributions that are suppressed by l/iV c 2 , where N c is the number of colors. When the emitted 
parton is a soft gluon, the splitting function should match the probability given by Feynman graphs in 
the limit p m +i — ► 0. Away from the soft and collinear limits, however, one can choose what functional 
form to use and one can debate the merits of different choices. 

This can be illustrated by the case of VlNClA [329], which represents a new development and is 
discussed in more detail later in this section. In a leading order shower in the leading color limit, the 
fundamental object that emits gluons is a color dipole, that is, two partons, say I and k, that are color- 
connected (i.e., adjoining on a color string). The basic idea here goes back to the Lund dipole [330], 
implemented in Ariadne [331]. (We shall henceforth refer to such showers as dipole-antenna showers, 
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in order to disambiguate them from what we shall call partitioned-dipole showers below.) The relevant 
Feynman graphs in the amplitude are those in which the gluon m + 1 is emitted from partem I and those in 
which it is emitted from parton k. In the squared amplitude, one has a contribution l-l, corresponding to 
the square of the graph for emission from I, a similar contribution k-k, and two interference contributions 
l-k and k-l. The approximation of keeping only the leading color contribution restricts us to the case that 
I and k are color-connected. In dipole-antenna showers, each dipole is treated as a unit, an antenna that 
radiates gluons, and the splitting functions can be chosen such as to match the perturbative result in all 
the relevant limits, i.e. gluon m + 1 collinear to I, collinear to k, or soft. There are two main differences 
between Vincia and Ariadne (and also a recent Sherpa implementation [332]). The first is that an 
explicit possibility to vary the shower ambiguities away from the singular regions is retained in VlNClA, 
and the second is that it combines the original dipole shower with the antenna subtraction formalism of 
Refs. [316, 318, 319] to match to fixed-order matrix elements. 

Another new development is what we can call the partitioned-dipole shower [333,334], which is 
discussed in more detail later in this section. Here one partitions the splitting function into two parts. One 
part contains the singularity corresponding to parton m + 1 being collinear with parton I and part of the 
soft singularity. The other part contains the singularity corresponding to parton m+1 being collinear with 
parton k, along with the remaining part of the soft singularity. Away from these singularities, one has a 
choice. A sensible choice (as suggested in Ref. [335]) is to take the splitting functions to be precisely 
those defined by the Catani-Seymour dipole subtraction scheme [175] that is widely used for next-to- 
leading order perturbative calculations. This has the advantage that it should be fairly straightforward 
to match these NLO calculations to a Catani-Seymour dipole shower. It has the disadvantage that the 
splitting of the emission probability from a dipole antenna into two parts is perhaps a bit artificial. There 
is more than one way to accomplish this splitting. 

The second part of the generator TL\{t) of shower evolution is the specification of the momentum 
mapping. In Eq. (140), there is a nominal integration over the momenta of all the partons after the 
splitting. However the matrix element of TL\(t) contains delta functions that, for given starting momenta 
{p}m, restrict the new momenta {p} m +i to lie on a three dimensional surface. This surface could be 
parametrized by splitting variables t,z,4>, as in Eq. (139). In the case of the timelike dipole-antenna 
showers in Ariadne, Vincia, and Sherpa, the momenta of all of the partons not part of the dipole 
remain the same before and after the parton emission. For the partons / and k that form the dipole, the 
momenta pi and p k plus three splitting variables are mapped reversibly to the momenta of three daughter 
partons, pi, p k , and p m +i after the splitting, with all of the parton momenta being on-shell. This mapping 
is symmetric under label interchange I *-t k. In the special case that p m +i is collinear with pi, we have 
pi = pi + pm+i and p k = p k . Similarly, if p m+ i is collinear with p k , we have p k = p k + p m+1 
and p\ = pi. In the soft limit, p rn +i = 0, we have pi = pi and p k = p k . Away from these limits, 
the mapping is necessarily not so trivial, leading to a further non-singular ambiguity which VlNClA 
attempts to explore. For the partioned-dipole shower, there is a similar but simpler mapping, this time 
not symmetric under I <-> k. The splitting function that includes the singularity for p m +i collinear with 
Pi, comes with a momentum mapping for which Pi = Pi + p m +i and p k = p k when p m +i is collinear 
with pi or soft. Away from these limits, the mapping takes some momentum from parton k in order 
to keep momentum conserved and all partons on shell. The choice here is to use the same momentum 
mapping as was defined by Catani and Seymour for the subtractions in next-to-leading order calculations. 
In the case of splittings involving an initial state splitting (spacelike showers), the momentum mappings 
are a little more complicated than sketched here. We should mention that it is also possible to take the 
momentum needed to keep all partons on-shell from all of the final state partons in what might be called 
a democratic way [328]. 

We hope that this comparative discussion may be useful as a guide to the more detailed presen- 
tations later in this section. We may also mention the published work [328] that was discussed at Les 
Houches but is not separately presented in this document. The idea here is to extend the idea of a lead- 
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ing order parton shower so that one does not average over spins or take just the leading color limit. In 
this case, there is an evolution equation similar to Eq. (140) but with spin indices and a more detailed 
specification of the color state. The solution of the evolution equation yields integrals that could, in 
principle, be computed numerically. However, an algorithm that is likely to be usefully convergent with 
finite computer resources is still under development [336]. 

17. TIME-LIKE SHOWERS BASED ON DIPOLE-ANTENNA RADIATION FUNCTIONS 30 

17.1 Introduction 

In this report we take the next step in the development of the VlNClA shower towards a full-fledged parton 
shower, embedded into the Pythia 8 generator [329,337]. Previously, we included only the gluonic 
time-like shower [329]. By including massless quarks we can start making comparisons at LEP energies 
and make quantitative studies for future linear colliders. As the VlNClA shower is a dipole-antenna 
shower, we can make direct comparisons with the dipole-antenna functions used in Ariadne [331]. 

We also make a phenomenological comparison with the Pythia 8 shower. For this purpose, we 
choose the evolution variable, the hadronization boundary and other parameters in VlNClA as close as 
possible to the default Pythia 8 settings. In this emulation mode we compare a few representative 
distributions, both infrared safe and infrared regulated observables, such as jet rates, thrust, and parton 
multiplicities for hadronic Z decays at yfs = mz- 

17.2 Dipole-antenna functions 

The most general form for a leading-log antenna function for massless parton splitting, ab — > arb, can 
be represented by a double Laurent series in the two branching invariants, 

^ oo 

a(y ar ,yrb; s) = - ^ C a> p y% r %. h , (144) 

S a,p=-i 

where 

s = s ab = s arb and y i:j = < 1 (145) 

are the invariant mass squared of the antenna and the scaled branching invariants, respectively. In prin- 
ciple, eq. (144) could also be multiplied by an overall phase space veto function, restricting the radiation 
to specific "sectors" of phase space, but we shall here use so-called "global" antenna functions which 
are summed together without such cuts. Note that we have here written the antenna function stripped of 
color factors, to emphasize that this part of the discussion is not limited to the leading-color limit. 

The coefficient of the most singular term, C_i _i, controls the strength of the double (soft) sin- 
gularity (the "double log" term) and the coefficients C_ij>o and Cj>o,-i govern the single (collinear) 
singularities ("single log" terms). These, in parton shower terminology collectively labeled "leading 
log" terms, are universal, whereas the polynomial coefficients Cj>Oj>o are arbitrary, corresponding to 
beyond-leading-log ambiguities in the shower or, equivalently, different NLO subtraction terms in the 
fixed-order expansion. 

We take the Gehrmann-de-Ridder-Glover ("GGG") antenna functions [319] as our starting point. 
The corresponding coefficients C a ^ for the the five antennae that occur in massless QCD at LL are 
collected in tab. 11. For reference, we also compare to the radiation functions [330, 338, 339] used 
in the Ariadne dipole shower [331], which are also the ones used in a recent study by the Sherpa 
group [332]. Note that the single log terms have a slight ambiguity when gluons are involved, arising 
from the arbitrary choice of how to decompose the radiation off the gluon into the two antennae it 
participates in. Nominally, the Ariadne single log coefficients therefore look different from the GGG 

30 Contributed by: R. Frederix, W.T. Giele, D.A. Kosower, RZ. Skands 
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Table 11: Laurent coefficients for massless LL QCD antennae (ab — > arb). The coefficients with at least one negative index 
are universal (apart from a re-parameterization ambiguity for gluons). For "GGG" (the defaults in VlNCIA), the finite terms 
correspond to the specific matrix elements considered in [319]. In particular, the qq antenna absorbs the tree-level Z — > qgq 
matrix element [320] and the gg antennae absorb the tree-level h° — > gg — > ggg and h° — > gg — > gqq matrix elements [322]. 
The qg antennae are derived from a neutralino decay process [321], 



ones. However, a re-parameterization of the total gluon radiation, which we label Ariadne2, reveals 
that the only real difference lies in the choice of finite terms. Interestingly, while all the Ariadne 
radiation functions are positive definite, the equivalent Ariadne2 one for gg — > ggg is not and hence 
could not be used as a basis for a shower Monte Carlo. 

In modern versions of Ariadne, gluon splitting to quarks has an additional pre-f actor 2/(1 + 
s-g/sg-), where c is the neighbor on the other side of the splitting gluon. This is based on comparisons 
to e + e~ — > qq'q'q matrix elements and implies that the smaller dipole takes the larger part of the g — > qq 
branching. Such effects are not included in VlNCIA at this point. 

Our convention for color factors is that they count color degrees of freedom. Their normalization 
should therefore be such that, in the large- Nc limit, they tend to Nc raised to the power of the number 
of new color lines created in the splitting. In particular, 

/V 2 1 

C F = ^ 



N C ~ 3 ' (146) 
C A = N c =3. 

For gluon splitting to quarks, the antenna shower explicitly sums over each flavor separately, hence the 
relevant antenna functions should be normalized to one flavor, Tr = 1. (We use the hatted symbols Cp 
and Tr to distinguish this normalization from the conventional parton-shower one in which Cp = 4/3 
and T R = 1/2.) 

The complete antenna functions, in the notation of [329, eqs. (2) and (11)], are then 
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where a s = a s (fiR) may depend on the branching kinematics. If so, we use a nominal a s = 1 for 
generating trial branchings, which are then accepted with probability a s ([iR) at the point when the full 
kinematics have been constructed (see below). The possibilities for [ir currently implemented in VlNClA 



are 



type 
type 1 
type 2 



K R 2 P± 

KrQe 
Kr 



(148) 



S ab ' 



where Kr is an arbitrary constant, p± is denned as in Ariadne with p\ = s ar s r b/s~^ [331], Qe is 
the evolution variable, and y^sTg is the invariant mass of the mother dipole-antenna. The default is a 
1-loop running five-flavor a s with /jr = p± (i.e., Type above, with Kr = |) and a s (mz) = 0.137 
(the default in Pythia 8, making comparisons simpler). Alternatively, both fixed and 2-loop running 
options are available as well [337]. For the pure shower, the dependence on the renormalization scheme 
of a s is beyond the required precision and hence we do not insist on an MS definition here. Indeed, the 
default value of a s (mz) in PYTHIA 8 is determined from tuning to LEP event shapes. Though beyond 
the scope of the present paper, we note that in the context of higher-order matching, one should settle 
on a specific scheme, and should then see the dependence on both the scheme and scale choices start to 
cancel as successive orders are included. 



17.3 Shower implementation 

Brief descriptions of the VlNClA switches and parameters are contained in the program's XML "man- 
ual", by default called Vincia . xml, which is included together with the code. This file also contains 
the default values and ranges for all adjustable parameters, which may subsequently be changed by the 
user in exactly the same way as for a standard Pythia 8 run [337]. 

The default antenna functions are contained in a separate XML file, Antennae-GGG . xml. An- 
tennae that are related by charge conjugation to the ones listed tab. 11 are obtained by simple swapping 
of invariants (e.g., gq antennae are obtained from the qg ones). Similarly, antenna functions that are 
permutations of the ones in tab. 11, such as gg — > qqg, are obtained by swapping. In view of the prob- 
abilistic nature of the shower, all antenna functions are checked for positivity during initialization. If 
negative regions are found, the constant term Co,o is increased to offset the difference and a warning is 
given, stating the new value of Co,o- 

We use the PYTHIA 8 event record [337], which includes Les Houches color tags [340,341] for 
representing color connections. At every point during the event evolution, leading-color antennae are 
spanned between all pairs of (non-decayed) partons for which the color tag of one matches the anti-color 
tag of the other. 

Shower generation proceeds largely as for the pure-gluon case described in [329], including the 
choice between two evolution variables 



VE 



type I (p_L-ordering) : y f = 9l = ± SarSrb = 4y ar y rb 



s s 2 



Q 2 ' (149) 

type II (dipole-mass-ordering) : yjj = — S = 2min(y ar , y r b) 



•s 



Note that we do not include an "angular-ordering" option. In conventional parton showers, which use 
collinear splitting functions, angular ordering gives a good approximation of the coherent dipole radiation 
patterns we here describe by the antenna functions A. Since dipole-antenna showers use A directly, 
coherence is thus independent of the choice of evolution variable to first order in this formulation (see, 
e.g., [330]). 

For the phase space map an optimal choice for the functional form of the "recoil angle" ijjaa 
(see [329,331]) away from the soft and collinear limit exists for qq antennae [342]. However, we have 
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not yet implemented this particular subtlety in the VlNClA code. The default choice for all antennae is 
thus currently the same as for the gg — > ggg splitting in Ariadne [331] 



RIADNE 



El + El 



(tt 



Jab) 



(150) 



with alternative choices listed in [329]. 

Trial branchings are generated by numerically solving for Atrial m tne equation R = A(y tria i), 
where R is a random number uniformly distributed between zero and one, and the trial Sudakov is [329, 
eq. (51)] 



trial J 



exp 



dye 

Atrial 
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1 ~ Var A(y, 
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with A an overestimate of the "true" antenna function such that 

A(y ar ,y rb ) = s arb A(^y ar , y rb , s arb , 1) > s arb A{y ar , y rb , s arb , 1) 



(151) 



(152) 



only depends on the rescaled invariants (for instance by using a fixed overestimate of a s = 1 here). Once 
the full kinematics are known (see below) the trial branching can be vetoed with probability 1 — A/ A, 
which by the veto algorithm changes the resulting distribution back to that of A, as desired. 

During program execution, cubic splines of A and A^ 1 are used for the actual trial generation. 
These splines are constructed on the fly, with the 2-dimensional phase space integrals in eq. (151) earned 
out either by 2-dimensional adaptive Gaussian quadrature (AGQ) on A directly or (substantially faster) 
by 1 -dimensional AGQ on the primitive function along a contour of fixed y ar , defined by 
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where a s is the overestimate of a s discussed earlier, Ci represents the color factors appearing in eq. (147), 
and the phase space limits depend on the choice of evolution variable, see below. During initializa- 
tion, the program checks for consistency between the analytic and numeric integrals and a warning is 
issued if the numerical precision test fails. 

The antenna with the largest trial scale is then selected for further inspection. A <p angle distributed 
uniformly in [0, 2ir] is generated, and a complementary phase space invariant, z, is chosen according to 
the probability distribution 
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where | J(yE, z)\ is the Jacobian arising from translating {y ar , y r b} to {yE, z} and z m i n (yE) is the small- 
est value z attains inside the physical phase space for a given yE- Depending on the type of evolution 
variable, as defined in eq. (149), we choose {yE, z}(y ar ,y r b) as 



type I : y E = Ay ar y rb , z = y rb 

=> |Jl| = l/(4z) , 2W,min(y£) = 

type II : y E = 2y ar , z = y rb for z < 1 

Ve = tyrb , z = y ar + (1 - 2y rb ) for z>l - \y E 



(155) 
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where, for type II, we have arranged the two separate branches y ar < y r b and y r b < y ar one after 
the other by a trivial parallel displacement in the z coordinate. Using the Laurent representation of the 
antenna functions, the analytical forms of I z become 
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where the I a is defined in eq. (153) and I* is the primitive along a direction of fixed y r b 
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17.4 Numerical results 

We now turn to a quantitative comparison between Pythia 8 and VlNClA for e + e~ — ► Z — > qq at 
= mz- We use a 1-loop running a s with a s (m^) = 0.137 (the default in Pythia 8), with a 5-flavor 
running matched to 4 and 3 flavors at the b and c thresholds, but to eliminate the question of explicit 
quark mass effects we only allow d and u quarks in the Z decay and subsequent shower evolution. 
The evolution is terminated at pj_had = 0.5 GeV, and we have switched off hadronization so as not to 
unintentionally obscure the differences between the partonic evolutions. Likewise, photon radiation is 
switched off in all cases, and in Pythia 8 we further switch off gluon polarization effects. For VlNClA, 
we use three different settings: transverse-momentum ordering with "GGG" antenna functions, dipole- 
mass ordering with "GGG" antenna functions, and transverse-momentum ordering with the "Ariadne" 
antenna functions. 

Fig. 31 shows the 3-, 4-, and 5-jet inclusive fractions as functions of the logarithm of Durham hj>, 
using the default Pythia 8 Durham clustering algorithm [337]. In Pythia 8, the 3-jet rate (the set of 
curves furthest to the right) is matched to the tree-level 3-parton matrix element, whereas the GGG and 
Ariadne antenna functions in VlNClA reproduce it by construction. The general agreement on the 3-jet 
rate is therefore a basic validation of the qq — > qgq antenna implementation. Higher-order effects appear 
to make the mass-ordered VlNClA slightly softer, which we tentatively conclude is due to this variable 
favoring soft wide-angle radiation over high-p^ collinear radiation (as illustrated by fig. 2 in [329]). 

Similarly, the 4-jet fractions (the middle set of curves in fig. 31) test the qg antennae in VlNClA, 
with the GGG showers here slightly higher and the Ariadne one slightly lower, in agreement with the 
differences in qg antenna finite terms, cf. tab. 11. This trend becomes more pronounced in the 5-jet 
fraction, since also the gg — > ggg function in Ariadne is softer than GGG. 

We may now study further distributions, as a representative example of which we take thrust, 
illustrated as 1 — T in the top row of fig. 32. The full distribution is shown to the left with a closeup 
of the region 1 — T < 0.1 to the right. The region 0.1 < 1 — T<|is dominated by well-separated 
three-jet configurations. In the tail, 1 — T > |, a matching to e + e~ — ► 4 jets would be required to 
improve the accuracy. In the region below 1 — T = 0.1, however, this would not help. These are three- 
jet configurations which are "nearly two-jet". Here, the type and size of the Sudakov suppression is 
essential, the first fixed order of which could be accessed by 1-loop matching, but since the fixed-order 
expansion is poorly convergent in this region anyway, the disagreement is more likely to be cured by a 
systematic inclusion of higher-logarithmic effects in the showers (either implicitly, by "clever choices" 
of evolution, renormalization, and kinematic variables in the LL shower, or explicitly, by a systematic 
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inclusion of NLL splittings). It should be noted, however, that hadronization and hadron decay effects 
are important in the region below 

1 - T ~ 1 - max(^) = mm( yij ) < ( A few GeV ) < 01 (160) 

m z 

where the x and y fractions pertain to 3-jet configurations. This complicates the separation of genuine 
non-trivial higher-log effects from non-perturbative effects when comparing to experimental data at cur- 
rently accessible collider energies. 

Finally, as illustration of an infrared sensitive quantity, in the bottom row of fig. 32 we plot the 
probability distribution of the number of partons produced at the shower termination for each of the 
four models. The total number of partons is shown to the left and the number of quarks (not counting 
anti-quarks) to the right. The definitions of p± in Pythia and in Vincia/Ariadne, respectively, are 
not exactly identical, but they have the same infrared limiting behavior [343], and hence a comparison 
of the number of resolved partons with a cutoff at p±had = 0.5 GeV should be meaningful. Since we 
have also chosen the same a s values etc., the basic agreement between the models in the lower left- 
hand plot in fig. 32 reconfirms that there are no large differences between the showers, even at the 
infrared sensitive level. Ariadne produces somewhat fewer partons, consistent with the Ariadne 
radiation functions being slightly softer. On the right-hand plot, however, it is interesting to note the 
first substantial difference between Pythia 8 and the Vincia showers. The Pythia shower produces 
significantly fewer quarks than any of the VlNClA showers, despite its being higher or comparable on 
the total number of partons (cf. the left-hand plot). A similar difference between parton and dipole- 
antenna showers was observed in an earlier Ariadne study [339], in which a comparison was made to 
the virtuality-ordering of traditional parton showers. It is interesting that we here observe the same trend 
when comparing to the Pythia 8 shower which is ordered in p±. Finally, we note that this difference will 
also have practical consequences; in the context of tuning of hadronization models, the VlNClA showers 
will presumably need a stronger suppression of non-perturbative strangeness production to make up for 
the larger perturbative production rate, as compared to Pythia 8. 
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Fig. 32: Top row: Thrust, 1 — T, Bottom row: Number of partons (left) and number of quarks (right) at shower termination, 
with 2 massless quark flavors. 



17.5 Conclusions 

We have presented the inclusion of massless quarks into the VINCIA shower algorithm, implemented as a 
plug-in to the PYTHIA 8 event generator. The dipole-antenna radiation functions are expressed as double 
Laurent series in the branching invariants, with user-specifiable coefficients. At the analytical level, 
we compare the coefficients of the "GGG" antenna functions [319] used by default in VINCIA to the 
Ariadne ones [331]. Modulo a re-parameterization of emissions from gluons, we find the double and 
single log coefficients to be identical, as expected. The finite terms, however, are generally somewhat 
smaller for the Ariadne functions. This represents a genuine shower ambiguity which can only be 
systematically addressed by matching to fixed-order matrix elements. 

At the phenomenological level, we have also compared to the hybrid parton-dipole shower in 
Pythia 8 [337] for e + e~ — > Z — ► qq at </s = mz- We find a good overall agreement, even at the level 
of an infrared sensitive quantity such as the final number of partons. For the number of quarks produced, 
however, Pythia 8 is markedly lower than any of the Vincia showers we have compared to here. 
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Fig. 33: Factorization-scale (/if) dependence of the total cross section for the W -boson production at LHC. An apparent 
fiF dependence of the inclusive W (W + jet) production cross section (open circles) is greatly reduced when we add the 
LLL-subtracted W + 1 jet production cross section. The summed cross sections are shown with filled circles. 
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18. LLL SUBTRACTION AND PS KINEMATICS 3 1 
18.1 Introduction 

We are developing NLO event generators for hadron collision interactions based on GRACE [344], using 
the Limited Leading-Log (LLL) subtraction technique [345] for the parton radiation matching. The 
matching technique is crucial since the contributions of an additional QCD parton radiation in NLO are 
also involved in the evolution of Parton Distribution Functions (PDFs) in a collinear approximation. A 
naive application of a PDF to NLO calculations results in an apparent double-counting. We avoid the 
double-counting by subtracting Leading-Log (LL) collinear contributions from the matrix element (ME) 
of radiative processes. The subtraction is stopped ("limited") at the factorization scale ((J,f) since PDFs 
do not involve any radiation harder than this energy scale. The LL contribution of the radiation is easy 
to calculate [346], though an appropriate care is necessary in the kinematical mapping to non-radiative 
processes [345]. The subtracted LL terms are formally moved to non-radiative processes and to be 
cancelled with divergences in virtual corrections. 

Figure 33 shows the sum of the total cross sections for inclusive W -boson production and LLL- 
subtracted W + 1 jet production evaluated for the LHC condition (proton-proton collisions at yfs = 14 
TeV). Here, "jet" denotes a gluon or a light quark in the final state. The cross sections are calculated 
using the tree-level MEs for W production and W + 1 jet production, respectively, convoluted with the 
CTEQ5L PDF [347]. Results are shown as a function of the factorization scale (/ip). We can see a strong 
fj, f dependence of the inclusive W production cross section (open circles) is greatly reduced by adding 
the LLL-subtracted radiative cross section. This shows a good matching between the ME and PDF; i.e., 
the LL contents in ME and PDF are nearly the same. 

31 Contributed by: S. Odaka 
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The virtual corrections are yet to be included in the results shown in Fig. 33. They can also be 
evaluated automatically in the framework of GRACE [346]. Divergent terms in these collections are to be 
cancelled with those moved from radiative processes. Remaining finite terms will alter the normalization 
of non-radiative processes, and will result in a substantial mismatch since there is no such correction 
in radiative processes. However, this mismatch is at the level of NNLO. It will be possible to restore 
the matching within the accuracy of NLO. The simplest way would be to change the normalization of 
LL components of the hard radiation remaining in radiative processes by the same amount as applied to 
non-radiative processes. This is actually a modification at the NNLO (a 2 ) level. 

So far we have discussed the matching in the integrated cross section. We have to achieve a good 
matching in differential cross sections, as well, in order to construct practical event generators. The QCD 
evolution evaluated in PDFs is simulated by means of a partem shower (PS) in event generators for hadron 
collisions. PS and PDF are based on the same factorization theory. However, since theoretical arguments 
are given only at the collinear limit, the theory gives us predictions only at the first-order approximation 
for the transverse behavior. It is necessary to introduce a model of 3-dimensional kinematics in order 
to construct a practical PS conserving the energy and momenta. The introduction of a suitable model is 
crucial for achieving a good matching in differential cross sections. We discuss about such models in the 
following sections. 



18.2 Initial-state PS kinematics 

We have constructed an initial-state Leading-Log (LL) PS program for the use in NLO event generation. 
The program is based on the simplest expression of the LL Sudakov form factor employing Q 2 as the 
evolving parameter, 



S(Qi,Ql) = «p 



« 2 dz ^m P(z) 



2 Q 2 J e 2vr 



(161) 



The details are described in our paper [345]. We stay in a naive LL implementation without introducing 
corrections partially incorporating higher order effects, such as the angular ordering, because we plan to 
extend our PS to a true Next-to-Leading-Log (NLL) approximation [348]. 

We first tested the kinematics model employed in the "old" PYHTIA-PS [349, 350], since the 
theoretical bases is nearly the same. We found this model gives a very soft transverse activity. It results 
in an apparent mismatch in the transverse momentum (px) distribution of W bosons, when we tried to 
merge the inclusive W production with the LLL-subtracted W + 1 jet production by applying this PS to 
both processes. The starting assumptions of the "old" PYHTIA-PS kinematics are that the z parameter 
of a branch is the ratio of squared cm energies after and before each branch instead of the fraction of 
light-cone momenta, and that the Q 2 is identical to the virtuality of the evolving partons. The first 
assumption requires the definition of a "target" partem; thus, it is model dependent. However, this z 
definition ensures a simple relation between squared cm energies of a hard interaction and the beam 
collision; Shard = ^l^-Sbeam where x\ and xi are given by the product of all z values in each beam. 
From a simple kinematical argument we found this model gives a relation, 

p 2 T = (1 - zfQ 2 , (162) 

for each branch at the soft limit [345]. 

On the other hand, ordinary arguments based on the massless approximation give a slightly differ- 
ent relation, 

p\ = (1 - z)Q 2 (163) 

at the soft limit. Apparently Eq. (162) gives a smaller p? value than Eq. (163) for a given set of Q 2 and 
z. The relation (163) must be better for the matching since external partons are nearly massless in ME 
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Fig. 34: Sum of the simulated pr spectra of W bosons for the inclusive W production and the LLL-subtracted W + 1 jet 
production at LHC. The new pr-prefixed PS described in the text is applied to both processes. Results are plotted for three 
different choices of /if: /-If/wiw = 0.5 (open circles), 1.0 (filled circles) and 1.5 (open squares). 



calculations. We have introduced a new kinematics model where px of each branch is given ("prefixed") 
by Eq. (163). We keep the definition of the z parameter. The momenta of evolving partons are calculated 
from this px value and the z value. Thus, the virtuality is not necessarily equal to the Q 2 of a branch. This 
new PS gives a harder W -boson pr spectrum than the "old" PYTHIA-PS in the inclusive W production 
simulation, showing a better matching to the LLL-subtracted W + 1 jet simulation. The sum of the two 
simulations gives a smooth px spectrum stable against a variation of the factorization scale {hf) [345]. 

After the submission of the paper [345] we tried another definition of the "prefixed" px, 

p T = (l-z- Q 2 /s)Q 2 . (164) 

The parameter s is the squared cm energy before the branch. This is the result of the massless approxi- 
mation of branching kinematics before taking the soft limit (Q 2 /s — ► 0). This definition is ugly in some 
sense since s is model dependent, but gives us a better matching than Eq. (163). We plot the summed px 
spectra of W -bosons for three different \ip values (jJLp/mw = 0.5, 1.0 and 1.5) in Fig. 34. We can see 
almost no variation of the spectrum except for a small difference around px = ray/ in this \i p range. 



18.3 Prospects for the final-state PS matching 

It is enough to consider the initial-state matching if we consentrate ourselves to NLO corrections for 
color singlet or heavy particle productions. However, once we go to NLO for those processes having a 
gluon or a light quark ("jet") in the final state, we also need to consider the matching in the final state. 

We plan to use a simple LL partem shower employing Q 2 as the evolving parameter also for the 
final state. We need to introduce an appropriate kinematics model to this PS, too. In the initial-state PS, 
models in which the definition of px precedes that of Q 2 give us a better matching as we have discussed in 
the previous section. This is because px is in principle an observable quantity while Q 2 is not physically 
well-defined for initial-state partons. Similar arguments should be done also for the final state. 

In the final state, the virtuality is in principle an observable as an invariant mass of particles even 
after the hadronization and decays. Therefore, it must be natural to identify Q 2 as the virtuality of the 
evolving partons. The transverse momentum (px) is also an observable in principle. Thus, z should 
be treated as an unphysical parameter. It should be used as a hidden parameter only to give px values 
according to the relation, 

p T = z(l - z)Q 2 . (165) 
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In this kinematics model, PS is a process to give final-state partons additional masses equal to the Q 2 of 
their first branches. The invariant mass of the hard interaction system should be unchanged even after 
the application of PS, since it is a very fundamental parameter for the evaluation of matrix elements. We 
also want to keep the production angles in the cm frame unchanged. These requirements can be fulfilled 
by introducing a common multiplication factor to the momenta of all final-state particles. 

We need to apply a proper mapping of non-radiative subsystem in a radiative event to an on-shell 
non-radiative event in the LL subtraction. A mapping using momenta of the branched parton and the 
target parton works well for the initial-state radiation [345]. The subsystem is boosted and rotated to its 
cm frame where the momenta of two incoming partons are aligned along the z axis. This is the process 
exactly reversing the kinematical rearrangement in our initial-state PS. 

The mapping should be done in the same concept also in the final state, exactly reversing the 
rearrangement in PS. It can be done as follows: pick up an arbitrary pair of final-state partons. If they 
can be considered as products of a PS branch, replace them with the parent parton having the invariant 
mass of the pair as its virtuality (Q 2 ). If not, skip this pair. Rescale the momenta of all particles in 
the cm frame with a common factor to make the replaced parent parton become on-shell. Evaluate the 
matrix element of the non-radiative process based on these rearranged momenta, multiply it with the LL 
radiation factor proportional to l/Q 2 , then we get an LL approximation of a final-state radiation. This 
procedure should be applied to all possible combinations if we have more than two partons in the final 
state. 

We expect that the LL contribution can be evaluated in such a systematic way, including the initial- 
state contributions, as well. All contributions should be summed to evaluate the total LL contribution. A 
program is under development based on these concepts. 

18.4 Conclusions 

We have achieved a good matching between PDF and matrix-element (ME) evaluations for the parton 
radiation in NLO QCD corrections, by using the Limited Leading-Log (LLL) subtraction technique. It 
has been demonstrated as a good stability of the W production cross section against a variation of the 
factorization scale (jif), where the total cross section is evaluated by the sum of the cross sections for 
inclusive W production and the LLL-subtracted W + 1 jet production. 

We have to achieve a good matching between the parton shower (PS) and ME, as well, in order 
to construct practical NLO event generators. The transverse activity of PS depends on the applied kine- 
matics model of parton branches. We have successfully built a suitable model for our Leading-Log (LL) 
initial-state PS, where p? is prefixed according to the relation in the massless approximation of branch- 
ing kinematics. The simulation employing this PS shows a good matching between the inclusive W 
production and the LLL-subtracted W + 1 jet production in the pt spectrum of W bosons. The spectrum 
is stable against the variation of /xf in a wide range. 

It is necessary to achieve a good PS -ME matching for the final-state radiation, as well, when we 
construct NLO event generators for those processes including "jet(s)" in the final state. A study is in 
progress for the final state based on the experience on the initial-state radiation. 
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19. A PARTON-SHOWER MODEL BASED ON CATANI-SEYMOUR DIPOLE FACTORISA- 
TION 32 

19.1 Introduction 

Parton-shower models form an indispensable building block of Monte Carlo event generators, such as 
Herwig [351], Pythia [352] and Sherpa [224], that aim at the realistic description of multi-particle final 
states as they are observed in high-energy collider experiments. By accounting for QCD bremsstrahlung 
processes, partem showers relate a small number of partons emerging from a hard interaction, denned 
at scale Qhard an d theoretically described through a fixed order calculation, to a larger set of partons 
at scales Q a <C Qhard- The parton-shower approach relies on the universal pattern of QCD emission 
processes once soft or collinear parton kinematics are considered. The soft and collinear phase-space 
regions are singular and obtain large corrections order by order in perturbation theory what makes an 
all-orders resummation of the associated kinematical logarithms essential. Most shower algorithms rely 
on collinear factorisation of QCD matrix elements and are accurate to the leading-logarithmic level. The 
Ariadne approach, however, is based around the soft limits [331]. 

The parton-shower approach being perturbative it cannot be extended to arbitrary small scales but has 
to be stopped at some infrared cut-off scale Q Q > Aqcd- Below that scale event generators model the 
transition of QCD partons into the experimentally observed hadrons through non-perturbative hadronisa- 
tion models. In fact, only through the incorporation of parton showers these hadronisation models can be 
made universal or independent of the underlying hard process. This, however, assumes that perturbative 
QCD between scales Qhard an d Qo is appropriately described by the parton-shower model used. 

In the past few years there have been lots of major improvements related to parton-shower Monte Carlos. 
This includes the incorporation of exact multi-leg tree-level matrix elements for the description of the 
first few hardest emissions from a given hard process, know as "matrix element parton shower merging", 
see e.g. [353,354], or the consistent matching of next-to-leading order calculations with parton showers, 
know as "Monte Carlo at NLO", see for instance [45,355]. In addition the available shower algorithms 
of Herwig and Pythia have been revised and improved [343, 356]. 

Only very recently new shower algorithms emerged that are based on formalisms used to construct sub- 
traction terms that allow for a numerical cancellation of infrared singularities in NLO QCD calcula- 
tions [329, 332-335, 357]. There exist now implementations of such shower algorithms for two com- 
monly used subtraction schemes, the antenna subtraction method [316] and the Catani-Seymour dipole 
formalism [175,358]. Besides incorporating the last knowledge on the infrared behaviour of QCD ma- 
trix elements, these models should largely facilitate the matching with NLO calculations earned out in 
the respective scheme. In this note we briefly report on the construction of a parton-shower algorithm 
relying on Catani-Seymour subtraction that has more extensively been presented in [333]. 

19.2 The shower model 

The Catani-Seymour formalism provides all the ingredients to construct a local approximation to the 
real-correction matrix element in any QCD NLO calculation. These subtraction terms, that can be con- 
structed in a process-independent way, possess exactly the same infrared divergences as the real-emission 
correction, such that the difference of the two is infrared finite and can safely be (numerically) integrated 
in four dimensions. In addition, the subtraction terms are chosen such, that they can be analytically inte- 
grated in d = 4 — 2e dimensions over the phase space of the produced soft or collinear parton that causes 
the divergences. The occurring 1/e 2 and 1/e poles exactly cancel the ones from the loop integration in 
the virtual part when adding the two pieces. Such, the Catani-Seymour method provides a way to con- 
struct a parton-level Monte Carlo program for a NLO calculations once the one-loop and real-emission 
corrections to the Born process are known. 

32 Contributed by: S. Schumann, F. Krauss 
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In the Catani-Seymour approach the additional soft or collinear partem is emitted from an emitter- 
spectator pair (called dipole). Considering both the emitter and the spectator to be either in the final 
or initial state, four configurations have to be considered, representing the singularities associated to 
emissions from the final or initial state. Labelling final-state particles by i,j and k and initial-state par- 
tons by a and b the real-emission matrix element can always be approximated by the sum over all the 
possible dipoles, 



i,j i k^i i 



(166) 



Hereby, T>ij^ describe splittings of a final-state parton ij into the pair i,j accompanied by a spectator 
k. Due to the presence of the spectator, four-momentum conservation and on-shell momenta can be 
accomplished locally for each individual splitting. The terms T>f- represent final-state splittings with an 
initial-state spectator, while and T> ai,b correspond to a splitting initial-state line accompanied by a 
final- and initial-state parton, respectively. The individual dipole terms are constructed from the Born 
matrix element by inserting colour- and spin-dependent operators that describe the actual splitting. For 
massless final-state emitters and final-state spectators, for instance, the dipole contributions read 

1 ~ T/^ T ■ ~ 

T^ij,k = ~~ ~z m\lj ■ ■ ■ , ij ■ ■ ■ , k, . . . =2 " ij ,k P-> • • • > *J • • • > • • • ) m ■ (167) 
2piPj Ttj 

The Tij and thereby denote the colour charge operators of the emitter and spectator, respectively, they 
lead to colour correlations in the full amplitude. The Vy^ are d-dimensional matrices in the emitter's 
spin space that induce spin correlations. 

For the construction of a parton-shower algorithm from the dipole formula Eq. (166) certain approxima- 
tions are needed that finally allow for an exponentiation of the splitting operators to derive the Sudakov 
form factors central for a shower implementation. In addition, the splitting kinematics, choices on scale 
settings and the actual shower-ordering parameter have to be fixed. 



19.2.1 Shower construction criteria 

The full colour correlations present in the |.A/f m+ i| 2 matrix element have to be discarded in the shower 
picture, instead the leading terms in 1/N C sue considered only 33 . In this approximation a colour flow 
can be assigned to each parton configuration. Motivated by considerations on the colour dynamics for 
soft emissions, we choose the emitter and spectator to be colour connected in the shower formalism. The 
colour-charge operators simplify to 

T k ' T ij 1 nfS1 

rp2 \f S ? eC ' V 1UC V 

ij J ^ 

with J\f^ ec = 1, 2 in case the emitter has one (SU(3) (anti-)triplet) or two (SU (3) octet) possible spec- 
tators. The four-dimensional dipole functions V are used as the shower splitting functions. Furthermore, 
we neglect spin correlations by using spin-averaged splitting functions (V) 34 . 

As shower evolution variable we choose the transverse momentum between the splitting products for 
branching final-state partons and the transverse momentum with respect to the beam for emissions from 
the initial state, collectively denoted by k^. This scale is also employed as the scale of the running 
coupling and the parton distributions, once initial-state partons are present. 



33 Although formally subleading, we consider splittings of the type g — > qq as well 
34 Some of the dipole functions can become negative in r 
interpretation. We choose to set them to zero in these cases. 



34 Some of the dipole functions can become negative in non-singular phase-space region, prohibiting a simple probabilistic 
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Based on the above approximations and choices Sudakov form factors corresponding to the different 
types of Catani-Seymour dipoles can be derived, that determine the probability for a certain branching 
not to occur for a given range of the evolution variable. The four generic cases are briefly reviewed in the 
following. For simplicity, here we consider massless partons only, the massive case is discussed in [333]. 



19.2.2 Final-state emitter — final-state spectator 

Consider the final-state splitting {ij,k} — > {i,j,k} with the four-momentum constraint pij + pk = 

Pi +Pj J rPk = Q and all momenta being on their mass-shell. The branching can be characterised by the 
Lorentz invariant variables 

PiPj - .. ~ PiPk , t ™ 

Vij,k = ; — , Zi = 1 - Zj = ■ . (169) 

PiPj + PiPk + PjPk PiPk + PjPk 

The factorised form of the fully differential (m + 1) -partem cross section that exactly reproduces the 
corresponding soft and collinear divergences of the real-emission process reads 

d<7 m+ i = da m ^2 ^2 dli V T~ \X ec ^ ~ yij,k)(V~ij,k{zi, yij,k)) ■ (170) 
The spin-averaged splitting kernels (Vijjt) for the branchings q — > qg, g — > gg and g —>■ qq read 

(v Wj ,fc(iy) = c F j _ - 2 (ni) 

(V gigj , k (zi, yij,k)) = 2C A { ,_-]_- + ^— — 2 + Zi (1 - %) \ , (172) 

( J- Zi + Ziyij k Zi + Hij^k ZiHij^ J 

(V mhk {zi)) = T R {1 - 2zi (1 - zi)} . (173) 

In terms of the splitting variables the transverse momentum between the splitting products i and j (our 
shower evolution variable) can then be written as 

k l = 2 PijPk Vij,k Zi (1 - Zi) , (174) 

and accordingly 

dVij,k _ dk^ 



Uij,k k 



(175) 



Setting the infrared shower cut-off equal to k^ and the upper limit to k^ max the z% integration is 
constrained to 



z T(kj_ imaX ' kj_ >0 ) — 2 I ^ ^ 



The kinematics of the splitting are fixed through 



\ 



k 2 

K _L,0 

k 2 

l.max 



(176) 



k 2 

Pi = ZiPij^ z — Pk + k±, (177) 

Zi 2pijp k 

Pj = (l-^ + ^_^_ A -fc ±J (178) 

Pk = 0--yij,k)Pk, (179) 
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with kj_ the spacelike transverse-momentum vector perpendicular to pij and p k and k±-k± = — \t 2 ± . The 
Sudakov form factor for having no final-state splitting with a final-state spectator between max and 



reads 



A FF(k^ ,k^_ 



exp 



EE \f si 



J\f spec 



V 



max -2+ 

2 



dk^ 



\ 



. (180) 



19.2.3 Final-state emitter - initial-state spectator 

In the presence of initial-state partons a final-state splitter may be colour connected to one of the incoming 
lines. We consider the splitting {ij, a} — > {i,j, a}, with pij — p a = Pi + Pj — Pa = Q ■ This time the 
branching is parameterised by the quantities 



PiPa + PjPa ~ PiPj 



J i],a 



1 — Zn 



PiPa 



PiPa + PjPa " PiPa + PjPa 

The relative transverse momentum of the new emerging final-state partons is given by 

1 — Xi 



ZpaPij 



Zi (1 - Zi) . 



(181) 



(182) 



The derived Sudakov form factor for this splitting type reads 

AFi(ki )max ,k[ )0 ) 
/ 

-EE 



cxp 



,max 

dk ± f Qs(kj) fajlla/Xjj^kl) 

- - - - n - \ v •l>ij,a)l 



\ 



Kj-spec 
ij a % 3 



k 2 



k 2 

"■±,0 



27T f a {Va,kl) 



(183) 



Here, i] a is the momentum fraction of the spectator parton a and f a (r] a , k 2 ) the coiTesponding hadronic 
PDF evaluated at some scale p? F = k^. The parton-distribution function f a {Va/xij,a, kj_) accounts for 
the new incoming momentum. The z,- L integration boundaries are given by Eq. (176) and the concrete 
splitting functions, (V^ (5j, Xij >a )), can be found in Ref. [333]. The branching kinematics are fixed to 



Pi 



Pj 



%i Pij ~r 



(1 - Zi)pij + — 



%i ^PijPa 
k 2 , 



Zi) 1p~ijPa 



Pa + k ± , 
Pa~k±, 



(184) 
(185) 



with k± perpendicular to both the emitter and the spectator momentum. The new spectator momentum 
is given by 



Pa 



-Pa 



(186) 



19.2.4 Initial-state emitter — final-state spectator 

Once a final-state line is colour connected to the initial state, besides the situation discussed in 
Sec. 19.2.3, the reversed case occurs as well. Namely, the initial-state line can split and emit a new 
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final-state partem while the spectator is in the final state. The momentum-conservation condition for such 
a branching {ai, k} — > {a, i, k} reads pf. — p a i = Pi + Pk — Pa = Q ■ The splitting variables are defined 
as 

PiPa + PkPa ~ PiPk 



x ik,a 



U: 



PiPa 



PiPa + PkPa PiPa + PkPa 

and the transverse-momentum squared of parton i with respect to the beam becomes 

i 2 o ~ ~ 1 — x ik,a /., \ 
k ± = 2p ai p k Ui{l - Ui) . 

x ik,a 

The Sudakov form factor associated with this splitting type reads 



(187) 



(188) 



A w (k] 



,max> *--L,0) 



exp 



EE \r sp 



Ki ec 

at k UL , 2 



max ^ 

2 



dk^ 



& x ik,a ^~ J{ x ik,ai u i'i i x ik,a, u i)) 



2tt 



with x_ = r] ai and x + = Q 2 /{Q 2 + 4k^ ) and 



?/ _ ,2 \ 1 u i ^ fa( 7 ]ai/ x ik,aikj_) 

•J { x ik,ai u i] — -, _ 9 7 7 O"^ 

1 ZUj X{k : a JaiyVaii ^-_}_) 



\ 
/ 

(189) 
(190) 



accounting for a possible flavour change of the incoming line through the backward-evolution step. The 
complete list of splitting kernels can again be found in Ref. [333]. The branching kinematics are given 
by 

Pa = —^—Pai, (191) 



Pi 

Pk 



x ik,a 

/1 \ x ik,a ~ 

(1 - Ui) = Pai + Uip k + fe_L , 



Ui 



x ik,a 
1 x ik,a ~ 



x ik,a 



Pai + (1 - - k± . 



(192) 
(193) 



19.2.5 Initial-state emitter - initial-state spectator 

The last case to be considered is the splitting of an initial-state line that is colour connected to the second 
incoming parton. The branching is parametrised through 



x i,ab 



PaPb ~ PiPa ~ PiPb 



V; 



PiPa 



PaPb PaPb 

such that the transverse-momentum squared of the new final-state parton becomes 

1 x i,ab Vi 



k ± = 2p ai p b Vi 

The Sudakov form factor of this configuration reads 
Anfk? 



x i,ab 



r 2 k 2 1 



exp 



EE 

j^spec 
ai b^ai a% 



k-! 



k 2 



(194) 



(195) 



r da nh J(x i>a ,, ki) (V ai ' h (x i>a ,)) 



2tt 



(196) 
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with 

J(x hab , Vl ;k ± ) = — , (197) 

and x_ = 7] a i and x+ = 2p a pb/(2p a pb + 4k^_ ). For the kinematics of the emission process it is 
convenient to keep the spectator momentum fixed and to align the new incoming parton a with the old 
incoming momentum according to p a = l/x^aft • Pai- The momentum of the newly emerged final-state 
parton i, is given by 

Pi = ~ — — —Pai + Vip b + k ± . (198) 

%i,ab 

Its transverse momentum has to be balanced by the entire set of final-state particles of the m-parton 
process (including all non-QCD particles). 



19.2.6 The algorithm 

Having at hand factorised expressions for all possible emission processes and corresponding Sudakov 
form factors a probabilistic shower algorithm of independent emissions can be formulated. The start 
seed forms a 2 — > 2 core event with fixed colour flow and a process dependent shower start scale k^_ max . 

1. The scale of the next emission is chosen according to the Sudakov form factors of all contributing 
emitter-spectator pairs. The dipole that yields the highest transverse momentum is picked to split. 

2. The value of the second splitting variable is chosen according to the splitting kernel. 

3. The splitting kinematics are determined, the new particle is inserted and the colour flow gets 
adapted. 

4. Start from step 1 as long as > and replace k^_ max by the transverse momentum of the 
last splitting. 

This yields a chain of subsequent emissions strictly ordered in transverse momenta. There is no formal 
subdivision of initial and final state evolution, instead, all dipoles are treated on equal footing. 



19.3 Comparison with experimental data 

The ultimate test of a theoretical model is a direct comparison with experimental measurements. Here 
we compare the newly developed and implemented parton-shower algorithm (called CS shower in the 
following) with some experimental data on hadron production in e + e~ annihilation, and Drell-Yan and 
jet production in pp collisions. Therefore the shower simulation has been supplemented with the string 
fragmentation routines of Pythia-6.2 [359] to account for hadronisation. 

We begin with some of the most precisely measured quantities, event-shape observables in e + e _ an- 
nihilation at the Z° pole. Fig. 19.3 contains a comparison for the normalised 1-Thrust (1 — T) and 
C-parameter (C) distributions with LEP1 Delphi data [360]. Both observables obtain large higher-order 
corrections for two-jet like events that appear as 1 — T « and C ps 0. In addition, there is a singularity 
in the C-parameter also in the region C « 0.75 that requires a resummation of large kinematical loga- 
rithms [324,361]. The CS shower yields a good agreement with the experimental data. Only very pencil 
like events, that are sensitive to hadronisation corrections, are overestimated in the Monte Carlo. We 
believe that this can be improved through a more detailed tuning of the hadronisation model parameters. 

In Fig. 19.3 we present the predictions of our model for the lepton-pair transverse-momentum distri- 
bution in Drell-Yan production and for the azimuthal decorrelation of inclusive dijet events in pp colli- 
sions. Both observables are nontrivial only if additional QCD radiation is produced and thereby test the 
emission pattern of the shower ansatz. We observe a good agreement with data for both observables in 
phase-space regions dominated by rather soft or collinear emissions but the agreement outside this range, 
i.e. large or small A^dijet* is also very satisfactory. 
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A6... distribution @ Tevatron Run II 




P T [GeV] ~ (rad) 



Fig. 36: The transverse-momentum distribution of e e Drell-Yan pairs compared with CDF Run I data [362] (left panel) and 
the Dijet Azimuthal decorrelation measured by D0 at Tevatron Run II [363] (right panel). 
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19.4 Conclusions 

We have presented a new parton-shower algorithm that uses fully factorised versions of the Catani- 
Seymour dipole functions to describe multi-parton production processes in a probabilistic manner. The 
model encodes exact four-momentum conservation on the level of each individual splitting due to the 
notion of splitting emitter-spectator pairs. Subsequent emissions are ordered in transverse momenta and 
the evolution of initial- and final-state partons is done in a unified way. Comparison with experimental 
data yields very encouraging results. In a next step we will combine this new shower approach with exact 
multi-leg tree-level matrix elements. Moreover, this model should facilitate a matching with exact NLO 
QCD calculations. 
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